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Traditionally high-aspect ratio triangular/tetrahedral meshes are avoided by CFD re- 
searchers in the vicinity of a solid wall, as it is known to reduce the accuracy of gradient 
computations in those regions. Although for certain complex geometries, the use of high- 
aspect ratio triangular/tetrahedral elements in the vicinity of a solid wall can be replaced 
by quadrilateral/prismatic elements, ability to use triangular/tetrahedral elements in such 
regions without any degradation in accuracy can be beneficial from a mesh generation 
point of view. The benefits also carry over to numerical frameworks such as the space- 
time conservation element and solution element (CESE), where simplex elements are the 
mandatory building blocks. With the requirement of the CESE method in mind, a rig- 
orous mathematical framework that clearly identifies the reason behind the difficulties in 
use of such high-aspect ratio simplex elements is formulated using two different approaches 
and presented here. Drawing insights from the analysis, a potential solution to avoid that 
pitfall is also provided as part of this work. Furthermore, through the use of numerical 
simulations of practical viscous problems involving high-Reynolds number flows, how the 
gradient evaluation procedures of the CESE framework can be effectively used to produce 
accurate and stable results on such high-aspect ratio simplex meshes is also showcased. 


I. Introduction 


In the multi-dimensional space-time conservation element and solution element’~° (CESE) method, tri- 
angles and tetrahedral mesh elements turn out to be the most natural building blocks for 2D and 3D spatial 
grids, respectively. As such, the CESE method is naturally compatible with the simplest 2D and 3D unstruc- 
tured grids and thus can be easily applied to solve problems with complex geometries. However, because 
(a) accurate solution of a high-Reynolds number flow field near a solid wall requires that the grid intervals 
along the direction normal to the wall be much finer than those in a direction parallel to the wall and, 
as such, the use of grid cells with extremely high aspect ratio (10° to 10°) may become mandatory, and 
(b) unlike quadrilateral/hexahedral grids, it is well-known that accuracy of gradient computations involving 
triangular/tetrahedral grids tends to deteriorate rapidly as cell aspect ratio increases, the use of triangular / 
tetrahedral grid cells near a solid wall has long been deemed impractical by CFD researchers.” 

In view of (a) the critical role played by triangular /tetrahedral grids in the CESE development, and (b) 
the importance of accurate resolution of high-Reynolds number flow field near a solid wall, a comprehensive 
and rigorous mathematical framework that clearly identifies the reasons behind the accuracy deterioration 
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as described above has been developed for the 2D case involving triangular cells. By avoiding the pitfalls 
identified by the 2D framework, and its 3D extension, it has been shown numerically that the accuracy 
deterioration phenomenon identified above can indeed be overcome completely. 


II. A Summary of Key Theoretical Results 


Here, we provide a summary of the key practical results obtained from the current mathematical devel- 
opment. The mathematically-inclined reader is referred to Section III for the more rigorous mathematical 
treatment and details. 


Figure 1. General notations of the vertices, internal angles, lengths of the sides and altitudes of a triangle. 


To proceed, let (i) AP, P:P3, depicted in Fig.1, be a triangle lying on the x — y plane with its area, 
A(AP, P2P3) >0 (2.1) 


and (ii) for each & = 1,2,3, let a, be the internal angle associated with the vertex, Py, I, be the length of 
the side facing the vertex P,, and hy be the length of the altitude originating from P,. Then, (i) 


1 
A= ale >0,k =1,2,3 2.2) 
(ii) 
I, > O and hy > 0,k = 1,2,3 2.3) 
and (iii) 
T > Q1,Q2,a3 > 0 anda; +a9+03 =7 2.4) 


Because of Eq.(2.3), the aspect ratio 7 of AP, P2P3 can be defined as 


a manta ln 13 
1) hi’ ha’ hs 


+>0 (2.5) 


i.e., 7 is the maximal value of the three positive ratios enclosed within the braces. In turn, by substituting 
a result of Eq.(2.2), ie., hy = 24, k= 1,2,3 into Eq.(2.5), one has 


lk 


ly i) l 1 ; ‘ é 
(2A/l1)’ (2A/l2)’ (Ain) = sqmaxt(li)”, (l2)"s (Is)"} > 0 (2.6) 


1 = max{ 


Let (ki, k2,k3) be a permutation of (1,2,3) such that In; > Ike > lp3. Then, with the aid of Eq.(2.3), we 
have 
lei = leo = les > 0 (2.7) 
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As an example, in the case where Ig > 1, > 13,k, = 2,k2 = 1, and k3 = 3. However in the case where 
1, = Ig > Is, one has either (i) ky = 1,k2 = 2, and k3 = 3, or (ii) ky = 2,kg = 1 and ks = 3. With the aid of 
Eqs.(2.2) and (2.7), Eq.(2.6) now implies that 


ea)? a? he 


= >0 2.8 
2A ler: Aer Pe oF) 


Note that Eqs.(2.2) and (2.7) imply that J,; and hy; respectively, are the lengths of the largest 
side and the shortest altitude of AP, P,P3. According to Eq.(2.8), the ratio of these two lengths 
is the aspect ratio 7 of AP; P2P3. 

To recast 7 as a function of the internal angles a,,a2, and a3, note that, by observing Fig.1 and using 
Eqs.(2.3) and (2.4), one has the following relations: 


hy, = lysinag = I3 sinag > 0 
ho = Is sina, = 1, sin az >0 (2.9) 


and hg = 1, sinag = lg sina, > 0 
Next, by using Eq.(2.9) and the definition of k 1, k2, and k3, one concludes that 
hry = le sin apg = 1x3 sinaz2 > 0 (2.10) 


In turn, by substituting Eq.(2.10) into Eq.(2.8), one has 


l l 
= Ee (2.11) 
[po SIN AK3 Ip3 SiN Ap 


Moreover, with the aid of a result of Eq.(2.4), ie., 
1> sina, >0,k = 1,2,3 (2.12) 


Eq.(2.3) and the law of the sines imply that there exists a constant coefficient 3 > 0 such that 


sin Qp1 sin Qp2 sin Ap3 


Le., 
lpi = Bsin Arp > 0, lho = Bsin ap > 0, and Ip3 = Bsin arn3 > 0 (2.14) 


As a result of Eqs.(2.7) , (2.12) and (2.14), we have 
1 > sinag > sinagzg > sinaz3 > 0 (2.15) 


Moreover, by substituting Eq.(2.14) into Eq.(2.11), one can cast 7 as a function of a%1, 2 and ayz3, i.€., 


(a2 (2.16) 
sin Qz2 SiN Ans 
As will be shown immediately, 7 can also be expressed as a function involving only az2 and ax3. 
With the aid of Eq.(2.4) and the definition of k,, k2,and k3, we have 
T > Ap, K2,%3 > 0 and apy + apg + OR3 = 7 (2.17) 
Thus, with the aid of some trigonometric identities, we have 
sin Qg1 = sin (7 — agg — OZ) = Sin (Ag2 + Ag3) = (Sin AK2)(COS AZ) + (Sin AK) (COS AK2) (2.18) 
Substituting Eq.(2.18) into Eq.(2.16) and using Eq.(2.17), we have 
N = cot age + cot args, O< aK2,QK3 <7 and Apo + AK3 <7 (2.19) 
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With the above preliminaries, the following proposition is proved in Appendix A. 
Proposition 1: Assuming Eq.(2.17), Eq.(2.15) implies that either 


> an > An. > ang >0 (case Eq.(2.17)) (2.20) 
or a T 
T > Ap > a and 5 > 1 — Ap > AK. > Ag > 0 (case Eq.(2.17)) (2.21) 


must be true. In turn, Eqs.(2.16) and (2.19)—(2.21) imply that: (i) 


def 2 

n= aE (2.22) 

(ii) 
7 = min if and only if ay = ag = a3 = 7 (2.23) 

(iii) 
li = 2.24 
Sie-s08 s 9 ( ) 

and (iv) 

n>>1 if and only if min{ay,a2,a3} = apg <1 (2.25) 


According to the above proposition, the aspect ratio 7 (i) attains its minimal value 2/\/3 if and 
only if AP, P,P; is equilateral, and (ii) attains a very large value if and only if the value of one 
of the internal angle of AP, P2P3; has a very small value. 

Next, let the “shape factor” - of AP, P2P3 be defined by 


1 © sin? ay + sin? ag + sin? a3 (2.26) 


It is shown in section III and Appendix B that (i) 


0<7< : (2.27) 
(ii) 
Y = Vmax 2 ; 1 = Q2 = 03 5  =k=l3 >0 (2.28) 
(iii) 
73 0t & mar{ay,a2,03} > 7 (2.29) 
and (iv) 
max{ay, 2,03} = Tex y=2 (2.30) 


2 


Hereafter, (i) the symbol “& ” is used to indicate that the statements given on its left and right sides are 
equivalent, while the symbol “ => ” is used to indicate that the statement given on its left side implies 
that given on the right side; (ii) as an example, y > 0* states that ~ approaches 0 from the range > 
0 while mar{a1,a2,a3} > a7 states that mar{a,,a2,a3} approaches 7 from the range < 7; (iii) as 
an example, the symbol {a1,a2,a3} denotes a set containing elements aj,a2, and a3, and, as a result, 
{a1, a2,a3} = {b1, b2, bs} indicates that the collection of elements aj, a2, and a3 is identical to that of the 
elements 6), b2, and bs; and (iv) as an example, (a1, a2, a3) = {a1, {a1, a2}, {a1, a2, a3}} so that (a1, a2, a3) = 
(b1, ba, bs) © ay = by, k = 1,2,3. 

According to Eqs.(2.27)—(2.30), (i) y attains its maximal value 9/4 if and only if AP, P2P3 is equilateral; 
(ii) y approaches its minimal value 0 if and only if the largest internal angle of AP, P:P3 approaches 7~ (and 
therefore its two other internal angles approach 0*), e.g., the highly obtuse triangle case depicted in Fig.2(a); 
and (iii) 7 ~ 2 and thus has a value relatively close to ymax if one of the internal angles of AP; P2P3 is close 
to the right angle even if one of its internal angles approaches 0 such as the case depicted in Fig.2(b). Note 
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P; 


OO, aa 


P; 
0 < a2, A3 «K 1 P; P; 
and 0<az, « 1,a,~> and 
0<n-a,<«1 a, ~m/2 


(a) (b) 


Figure 2. Representations of two kinds of high-aspect ratio cells: (a) highly obtuse, and (b) nearly right-angle triangle. 


that Ig/l3 = sinag/sina3 © ag/a3 if 0 < a2,a3 < 1. Thus, depending on the ratio a2/az, the ratio of 
the lengths of the two sides P,P; and P,P; can take any value > 0 even for the highly obtuse triangle case 
depicted in Fig.2(a). 

Even though the two triangles depicted in Fig.2 have shape factors of vastly different values — one of them 
approaches the minimal value while another is relatively close to the maximal value of 9/4, both triangles 
have large aspect ratio. This is due to the fact that both meet the condition that min{a,,a2,a3} « 1 
which, according to Eq.(2.25) is the necessary and sufficient condition for 7 > 1. 

In Section III, the relation between accuracy deterioration of gradient computation and the shape factor 
7 is established using a procedure briefly described below. First, for each k = 1,2,3, and at any given time 
level, let (i) 6, be the numerical value of a scalar function ¢ assigned to the vertex P, of a given triangle 
AP, P2P3 lying on the x — y plane, and (ii) (az, yz) be the coordinates of the point P,. Second, let Vd, 
the gradient vector of ¢ associated with AP, P:P3, be evaluated in terms of ¢,, 2%, yx, Kk = 1,2,3, using 
the standard linear interpolation procedure defined by Eqs. (3.2), (3.5), (3.6) and (3.21). Third, for each 
k = 1,2,3, let Ad, denote the numerical error of ¢, at the given time level introduced through a time- 
marching scheme while it is assumed that the coordinates (xx, yx), k = 1,2,3, be specified values with no 
numerical errors. Fourth, with the above assumptions it is shown in Eqs.(3.24), (3.35)-(3.37) and (3.40) 


that the numerical errors €,,€2 and e€3 of the directional derivatives along the side directions P)P3, P3P\, 
and PP of AP; P2P3, respectively, are dependent only on [x and Adz, k = 1, 2,3, while, according to Eqs. 
(3.26), (3.40)—(3.43), (3.46), (3.56), (3.57) and (3.66)—(3.69), the numerical errors Av, and Av, of 0¢/0x 
and 0¢/0y, respectively, are functions of Adz, x, and yz, k = 1,2,3 only. Fifth, the error norm for €1,€: and 
€3 is defined as \/(e? + 6} + €4)/3 while that for Av, and Avy are defined as \/[(Av,)? + (Av,)?]/2. Sixth, 
as shown in Eq.(3.208), the error amplification factor R for gradient computation in turn is defined as the 
ratio between the error norm for Av, and Av, and that for €;,€2 and es. Seventh, according to Eq.(3.209) 
and other associated equations given in Section III, it turns out that the factor R is a function of the internal 
angles a 1,@2 and a3, and the errors €1,€2 and €3 only. In fact, according to Eqs.(3.167)—(3.170), the 2 x 2 
real symmetric matrix H, which appears in Eq.(3.209) is a function of a1,a2 and a3 only, while, according 
to Eqs.(3.54), (3.108), (3.112), (3.115) and (3.118), the 2 x 1 real column matrix ~, which also appears in 
Eq.(3.209) is a function of a1, a2,a3,€2 and €3 (note: €1,€2,€3,Q1,Q@2 and ag are linked by Eq.(3.89), as 
such, given a set of a1,@2 and a3 , only two of €),€2 and €3 are independent parameters). Finally by using 
Eq.(3.209) and the Rayleigh-Ritz theorem(Ref.[8], p.43), one concludes that for, any combination of €1, €2 


and €3, (i) 
Vey) 5 RS 4/er(9),0<7 5974 (2.31) 
where, 
o+(7) st [3 +2 (074) - 7] 0<7<9/4 (2.32) 


are the two eigenvalues of H,, and (ii) each of the lower and upper bounds ,/a_(7) and ,/o4(y) in Eq.(2.31) 
can be attained by R by some special combinations of €1,€2 and €3. As such \/o_(7) and \/o+(7) are the 
greatest lower bound and the least upper bound of R, respectively. 

According to Eq.(2.32), both eigenvalues of H; are functions of the shape factor y only, even 
though H, itself is a complicated matrix function of a), a2 and a3 defined by Eqs.(3.167)—(3.170). 
Also, by using Eqs.(2.27) and (2.32), it is shown in Section II that, as the value of 7 decreases from its 
maximal value 9/4 to its minimal limit 0*, (i) the value of o;(7) increases monotonically from 
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1 to +00; and (ii) the value of o_(7) decreases monotonically from 1 toward the limit value of 
1/2. Thus, with the aid of Eqs.(2.28) and (2.29), one has (i) 


o(y=H=leo(y=l1eM=alh {the 2x 2 identity matrix 
9 es (2.33) 
Y 7 Qa, = a2 = a3 3 ly =lg =13 > 0 
(ii) 
1 ‘ 9 
5 <oa_(y)<1l<o(y), if0<7< ri (2.34) 
(iii) 
1 
a_(y) > Gy. & 04(7) 3 +00 6 ¥ 3 07 & mar{ay, 02,03} 3 1 (2.35) 


and (iv) 
3 
max{ay, 2,03} = 5 S7y=25 Voy) = ie e 1.225 (2.36) 


With the aid of Eqs.(2.33)—(2.36), Eq.(2.31) implies that, (i) the error amplification R for gradient 
computation has the constant value 1 for any combination of €,,¢€2 and ¢3, if and only if AP, PjP3 is 
equilateral - as such the most accurate gradient computation occurs in the case where the triangular grid 
is built from equilateral triangles, (ii) for the case in which the grid is built from triangles similar to that 
depicted in Fig.2(b), the least upper bound of R is only slightly greater than 1 - as such accuracy deterioration 
of gradient computation is rather slight for this case, and (iii) because \/o4(7) approaches +00 as y > OT, 
accuracy deterioration of gradient computation progressively becomes worse when triangles used in grid 
construction become more obtuse and it would reach an unacceptable level if these triangles become highly 
obtuse, similar to the triangle depicted in Fig.2(a). Fortunately, the triangular grids suitable for two- 
dimensional viscous simulations near a solid-wall boundary can always be constructed using triangles similar 
to that depicted in Fig.2(b). As will be shown in Section IV, for grids containing triangles similar to that 
depicted in Fig.2(b), the gradient evaluation procedures within the CESE method can be effectively used to 
improve numerical stability while maintaining numerical accuracy. 

Note that, as will be shown in a future paper, the mathematical development presented here for triangular 
grid can be extended in a straightforward manner to a tetrahedral-grid case, albeit the algebra becomes 
more complicated. In particular, the square of the error amplification factor R for a tetrahedral-grid case 
can still be cast into the form given on the extreme right side of Eq.(3.209) with the understanding that 
H, is a 3 x 3 real symmetric and positive-definite matrix [Ref.9, p.250] while ~ is a 3 x 1 real column 
matrix. In fact, it can be shown that, for a regular tetrahedron (i.e., a tetrahedron with all 
its four faces being equilateral triangles), H, is reduced to the 3 x 3 identity matrix, and as 
such the factor R = 1 for all possible combinations of the numerical errors associated with 
the directional derivatives evaluated along all six edge-directions of the tetrahedron. In other 
words, no accuracy deterioration of gradient computation occurs for a grid built from regular 
tetrahedrons. 

Moreoever, it can be shown that, for a tetrahedron in which three right internal angles share 
a common vertex (trirectangular tetrahedron), the least upper bound of R is 2 ~ 1.414, a 
number slightly larger than 1. As such accuracy deterioration of gradient computation will be 
mild for a grid constructed from such tetrahedrons. 


Ill. A Mathematical Framework for Identifying the Cause and Cure of 
Accuracy Deterioration Associated with Triangular Grids with High 
Apsect Ratio 


As a preliminary consider Fig.3(a). Here points P,, P2, and P3 are the vertices of a triangle lying on the 
x—y plane with their coordinates, relative to the Cartesian x — y coordinate system shown in Fig.3(b), being 
(21,41), (2, y2) , and (3, y3) respectively. Let (i) 


A(AP, P> P3) def the area of the triangle AP; P,P; > 0 (3.1) 
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90° counterclockwise 
rotation 


P, 


(1 ¥) 


XxX 
(a) (b) 


Figure 3. General representation of (a) the triangle, along with its vertices, coordinates and internal angles; and (b) 
Cartesian coordinate system. 


(ii) the determinant of a square matrix M be denoted by |M|_ , (iii) 


T2—-2%1 Y2- VY 
3-21 Y3-Y1 


= 


= (x2 — £1) (ys — yi) — (@3 — £1) (ye — y1) (3.2) 


and (iv) the rotation from the positive direction of the x—axis to that of y—axis is in the counterclockwise 
direction with an angle of 90°.Then it can be shown that: (i) 


for the case in which the boundary of AP, P:P3 traced in the P, + P,; — P3 — P, direction forms a 
counterclockwise loop, and (ii) 
6 =2-A(AP, P2P3) <0 (3.4) 


for the case in which the boundary of AP, P:P3 traced in the P, + P,; — P3 — P, direction forms a 
clockwise loop. 


P 


x 
Figure 4. Representation of the x — y — ¢ space. 


Next, for each k = 1, 2,3, let (i) ¢, be a scalar value assigned to the spatial point P,, and (ii) P, denote 
a point in the x — y — @ space depicted in Fig.4 with the coordinates (Xk, Yk, bk). In the following, first we 
will show that points Pi. Py, and Ps lie on the plane in the x — y — ¢ space represented by 


g=ax+byt+e (3.5) 


where, 
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TY 1 


g2-%1 yo" tg—-21 go—- hy tT Yo 2 
$3 —¢@ = t3-%, $3-—¢ 
qe eed Zi, pot He ‘I and oa aS vss (3.6) 
) 0) é6 
Proof: By definition, a plane in the x — y — ¢ space is represented by 
dz + doy + d3¢ + d4 =0 (3.7) 
where dj, dz, d3, and dy, are real constants with 

(di, dz, ds) 4 (0,0, 0) (3.8) 


i.e., at least one of dj,d2, and dg is non-zero. 
Next we will show that P,, P2, and P3 can only lie on a plane represented by Eq.(3.7) with dg 4 0. To 
prove by contradiction, let P;, P2, and P3 lie on a plane represented by Eq.(3.7) with dz = 0. Then 


dix, Tr doy. =I d4 =0 (3.9) 

dy x2 doy2 d4 =0 (3.10) 
and 

d\x3 d2y3 d4 =0 (3.11) 


By subtracting Eq.(3.9) from Eqs.(3.10) and (3.11), respectively, one has 


(x2 — 21)d, + (y2 — yi)d2 = 0 


3.12 
and (a3 — %1)d, + (y3 — yi)d2 = 0 : 
On the other hand, according to Eqs.(3.2)—(3.4), we have 
Sees aH (3.13) 
%3—%1 Y3— V1 


In turn, according to elementary algebra, Eqs.(3.12) and (3.13) = d, = dz = 0. Because Eq.(3.7) does not 
represent a plane in the x — y — @ space if dj = dz = d3 = 0, points P,, Py, and P3 can only lie in a plane 
represented by Eq. (3.7) with d3 4 0. 

Let d3 4 0. Then Eq.(3.7) © (i-e., is true if and only if) 


g=azrt+by+c (3.14) 


where d d d 
d@ 24 yy S_2 andd S-4 (3.15) 
ds 3 ds 
Note that, according to the above derivation of Eq.(3.14), the plane intersecting any three points Py, Py, 
and P3 with A(AP, P2P3) 4 0 (ie., Pi, Po, and Ps are not collinear points on the x — y plane) must be 
represented by Eq.(3.14) with the coefficients a’, b’, and c’ satisfying the conditions 


ar+byt+e = $4 
arg + b'yo +c = d2 (3.16) 
and a’x,+b'y3+¢ = ¢3 


Because (i) Eq.(3.16) represents a system of three linear equations for the coefficients a’, b’, and c’, and (ii) 
with the aid of Eq.(3.13), we have 


2, yw il Ly Yi 1 ieee Aid 
= _ (2-2, Yyeo-yi 
Z2 yo lj=|to-a21 yo-y. OJ = =d#40 (3.17) 
T3— 21 Y3— Yl 
zz y3 1 w3—-2, y—-yr O 
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according to elementary theory of linear system of equations, the coefficients a’, b’, and c’ can be uniquely 
determined in terms of xz, yx, and dp, k = 1,2,3. In fact, with the aid of Eqs.(3.6) and (3.17), we have (i) 


gi yi i P1 Yi 1 
2 yo 1 d2-¢1 yo-y O d2-o1 Yyo-y 
$3 y3 1 $3—¢1 ya-y O $3—¢1 Y3-Y1 
fies = = = 3.18 
. 4 5 5 s ue) 
(ii) 
xy gy 1 Beal or 1 
x2 d2 1 t2—-X, d2—-h, O 2-2, b2-h1 
rz $3 1 3-21 $3—¢, 0 3-21 63-1 
f= = - =b 3.19 
6 6 6 (3-19) 
and (iii) 
TZ YM 1 
r2 Yy2 2 
a X3 . 03 2% (3.20) 
Thus the plane defined by Eqs.(3.5) and (3.6) is the same plane defined by Eqs.(3.14) and (3.18)—(3.20). 
QED. 


Hereafter let the plane defined by Eq.(3.5) and (3.6) be denoted by [. Then, on I’, Eq.(3.5) implies that 
the gradient vector 
=, e O pas O > ay => 
Vee aes + a = aé; + bey (on T) (3.21) 
where é€;, and é, are the unit vectors in the c— and y— directions respectively. As such, on I’, (a) Vo isa 
constant vector; and (b) both 0¢/0x and 0¢/0y are constant coefficients. 
Next, for each k = 1, 2,3, let = 
Ph = Grex + Yrey, K=1,2,3 (3.22) 


ie., P, is the position vector of point P, on the x — y plane. Also, for each (k1,k2) with ky 4 kg and 
ky, ko = 1,2,3 ; let 


> de > > 2 “% 
Py, Pee = Pr, — Pry = (thy — ©, ES + (Yeo — Yer )G, er # ky and ky, ko =1,2,3 (3.23) 
and 
e > 
Sky ,ko ae Pry Pho = VJ (Lk =— Li)? + (Yio — Yk); ky x ko and ky, ko — 1,2,3 (3.24) 


—_— 
(Note here k; and kz are defined independent of Eq.(2.7)). According to Eq.(2.23), Py1 Pro is the displace- 
ment vector joining points P;,, and P,, and pointing towards P,, from Py, . Also by combining Eqs.(3.1) 
and (3.24), one concludes that 


Sky,ko = Sko,k1 > 0, ky # ko and ky, ke = 1, 2,3 (3.25) 
As a result, one can define 
as 
ef Pe, Pi e = 
Chg eK and «= fie ky Pho = Phe ky # kg and ky, ko = 1,2,3 (3.26) 
Sky ,ke Sky ,ke 


Moreover, Eqs.(3.23)—(3.26) = (a) each €;, 4, is an unit vector, i-e., 


|Ek, ke | = 1, ky x ko and ky, ke => 1,2,3 (3.27) 
and (b) 
Cho ,ky = —Cky,ko AN pk ,k, = Pky ko» ki A kg and ky, ko = 1,2,3 (3.28) 
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As an example, consider the plane [ which intersects points P,, Po, and P3. Let Vo be the constant 
gradient vector on P. Then Eqs.(3.21), and (3.23) > Vode Px, Pex = bk, — bk, Which, with the aid of Eq. 
(3.26) in turn > 


Vb © Eby ky = Why kgs ki # ky and ki,kg =1,2,3 (on plane L) (3.29) 
In other words, for each pair of k; and ky with ki # ko and ky, ko = 1,2,3, x, ,n, is the directional 
derivative of ¢ along the é;,,,,, direction. The value of p;,,,, can also be considered as the 


“slope” of P;, P,. in the «— y — ¢ space. 
Moreover, Eqs.(3.25), (3.26) and (3.28) imply 


—Ha1 = 1,2 = oi (3.30) 
$1,2 
— 13,2 = 12,3 = eae (3.31) 
$2.3 
and 
— 11,3 = 23,1 = ges (3.32) 
§3,1 
In turn, Eqs.(3.30)—(3.32) imply 
81,2 ° fl1,2 + $2.3 ° 12,3 + $3,1° U3,1 = 0 (3.33) 
and 
Haji + Hi2 = 13,2 + H2,3 = 11,3 + 3,1 = 0 (3.34) 


Because the line segments P; P2, P2P3 and P3P, are the three sides of AP; P2P3 , given the scalar values 
$1, G2 and os, f1,2, 12,1, 3,2; 2,3, 11,3 and 431 represent all possible directional derivatives of @ which can 
be evaluated along the spatial directions aligned with the three sides of AP, P2P3. Because the six directional 
derivatives are linked by four independent conditions given in Eqs.(3.33) and (3.34), any one of them can 
be determined in terms of any two independent directional derivatives among them. In fact, any two of 
them can be chosen as the independent directional derivatives, as long as they are evaluated 
along two different sides of AP,P:P3; and thus are not linked by one of the three conditions 
given in Eq.(3.34). 

For each k = 1, 2,3, let Ad, denote the numerical error of ¢, at some given time level introduced through 
a time-marching procedure. Moreover, (i) let the coordinates (x,, yz) of point Py,k = 1,2,3. be given 
fixed values that do not vary during the time-marching procedure, i.e., the numerical errors (Ax;,, Ay;,) of 
(tk, Yk), k = 1, 2,3, are assumed to be zero, and (ii) for any pair of ky and ky with ky 4 ko and ky, ko = 1, 2,3, 
let Apiy,,~. denote the error of 4x, ,~. introduced as a result of the errors Ady, and Ady, of o%, and dx, 
respectively. Then Eqs.(3.30)—(3.34) => 


Adz — A 
—Apiat = Apia = Bere: (3.35) 
$1,2 
Ads — A 
—Aps3,2 = Apie,3 = as (3.36) 
$2,3 
Ad, —A 
—Apii.3 — Aus,1 _ Agi Ags (3.37) 
$3,1 
81,2 ° Apo + 82,3 - Auo3 + 83,1 -Apu31 = 0 (3.38) 
and 
Ape, + Api = Apsz.2 + Ape,3 = Api3 + Aus. =0 (3.39) 
10 of 38 


American Institute of Aeronautics and Astronautics 


Because the six numerical errors Aug, .4.,k1 #4 ko and ky, kg = 1,2,3 are linked by four independent condi- 
tions given in Eqs.(3.38) and (3.39), any one of them can be determined in terms of any two independent 
errors among themselves. Any two of these errors can be chosen as independent errors as long as they 
are not linked by one of the three conditions given in Eq.(3.39), i.ec., the numerical errors of any two 
directional derivatives evaluated along two different sides of AP, P,P; can be chosen as the two 
independent numerical errors. 

To simplify the following developments, let 


€y dc Apes; €2 pie Apz,1 and €3 = Api2 (3.40) 


and 


~ def _, ~ def , ~ def , 
€1 = €2,3, €2 = €31 and e€3 = €1,2 (3.41) 


Then, with the aid of Eqs.(3.28) and (3.39), Eqs.(3.40) and (3.41) > 


Ausz,2 = —€1, A“i3 =—€2 and Ape1 =—€3 (3.42) 


and 


€3.2 = —é1, €13 = —€5 and €2,1 => —€3 (3.43) 


At this juncture, note that the integers 1, 2, 3, appear in Eq.(3.40) as component-identification indices. 
Moreover under the index permutation, 


(1, 2,3) > (2,3,1) (3.44) 
(i.e., the integers 1, 2, and 3, respectively are replaced by 2, 3, and 1), Eq.(3.40) will turn into the form 


€92 ple Aps.1; €3 ios Apii2 and €, as Api2.3 (3.40a) 


i.e., with the exception of them being presented in different order, the definitions given in Eq.(3.40)a are 
identical to those given in Eq.(3.40). Hereafter, an equation or a set of equations that possesses such a 
property is said to be invariant under the index permutation of Eq.(3.44). In fact, each similar definition set 
to be given in Eqs.(3.109)—(3.111) is also invariant under the same index permutation. 

Note that, because the numerical errors (Ax,, Ay,) of (x, yx),& = 1,2,3 are assumed to be zero in 
a time-marching procedure, according to Eqs.(3.22)—(3.26), for each pair of ky and ke with ki 4 kg and 
ky, kg = 1,2,3, the unit vector €,,, is a constant vector with no numerical error during a time-marching 
procedure. As such with the aid of Eqs.(3.5), (3.21) and (3.40)—(3.43), and the definitions: (i) 


Vy, = “0 =a and vy e a =b (onT) (3.45) 
and (ii) 
Av ees Avzé; + Avyéy (on T) (3.46) 


where Av, and Av,, respectively, denote the numerical errors of vy, and vy. One concludes from Eq.(3.29) 
that 


Ave & = €x,k =1,2,3 (onT) (3.47) 
Also, because of the assumption Eq.(3.1) , (i) Eqs.(3.26) and (3.41) > 


Ck, x tho if ky x ko and ky, ko = nl 2, 3 (3.48) 


and (ii) the internal angles aj, a2 and a3 of AP; P2P3 depicted in Fig.3(a) must satisfy the condition (Note: 
the symbol “ €” = “belongs to” ): 


(a1, Q2, a3) E Da (3.49) 


where 
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Ds a {(a1, a2, a3)| 0 < a1,02,a3 < 7 and a, + ag +a3 = 7} (3.50) 


Moreover, with the aid of Eq.(3.46), we have 


A> Sf VA Av = \/(Avz)? + (Avy)? (3.51) 


and thus 


(|ax)) = (Av,)? + (Avy)? (3.52) 


With the above preliminary, we will prove the following theorem: 
Theorem 1: Let (i) a1,a@2, and ag be the internal angles of AP, P)P3 , shown in Fig.3(a) and thus 
(Q1, A2, a3) € Da, (ii) 


aE 
S(a) & — ( : eal O<a<nr (3.53) 
sin* a \ cosa ‘I 
and (iii) 
12(2).0% (2) ane (2) o 
€3 €1 €2 
Then 
2 
(|Ax]) = (é,)'S(ax)éx, for each k = 1,2,3 (3.55) 


where for each k = 1, 2,3, (é,) is the transpose of é,. 


» 


> 


Ck 


0, k=1,2,3 
Xx 


Figure 5. Representation of a unit vector in x — y space. 


Proof: First note that, for each k = 1, 2,3, (i) & is the unit vector defined by Eqs.(3.26) and (3.41), and 
(ii) the angle 6; , shown in Fig.5, is uniquely defined by &, through the conditions: 


Ex = (cos Ox )Ez + (sin Ox )éy, k=1,2,3 (3.56) 
and 
t>0,>—7, k=1,2,3 (3.57) 
Also it follows from Eq.(3.56) that 


Ex, © Cx, = (COs Og, )(cos Ox.) + (sin Og, )(sin 0%, ) = cos(Ox, — Or, ), ki, ke = 1,2,3 (3.58) 
Because €,,k = 1,2,3, are unit vectors, Eqs.(3.48) and (3.58) => 


(Ex, @ Ek) = | cos(Ox, — On, )| < Lif ki A kg and ky, ke = 1,2,3 (3.59) 
which, in turn, > 
sin(O,, — Or,) AO if ky A ko and ky, kg =1,2,3 (3.60) 
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Note that, with the aid of Eq.(3.59), it can be shown that Eq.(3.48) <= Eq.(3.60). Moreover, because 


Eq.(3.57) > 
27 > Ok, — Ok, > -2r, if ky # kp and ky, ke = 1253 


As such, assuming Eq.(3.57), Eq.(3.48) <= Eq.(3.60) = 


Big — On, Fm and Op, — On, #0 if ky A ko and ky, ko =1,2,3 


Substituting Eqs.(3.46) and Eq.(3.56) into Eq.(3.47), one has 


(cos 6, )Av, + (sin@,)Av, = , for any k = 1,2,3 


On the other hand, given any pair of k, and ky with ky 4 kp and k1, ky = 1, 2,3, in turn Eq.(3.63) > 


(cos O,, )Avy + (sin Oy, Avy = €x, 


ky # ko and ky, ko = 1,2,3 
(cos 0%, Ave + tama 1# ke and ky, ke Ds 


Because, with the aid of Eq.(3.60) one has 


GOS ie SUE Salata Snape Nie (ode Oe ahi a= Sin Oa) 0 


cos 6, sin Oz, 


if ky x ko and ky, ke = 1,2,3 
Eq.(3.64) can be inverted to yield the following result 


6 = M(ky,ko)é(ki,k2) if ky A ke and k1, ko = 1,2,3 


where (i) 
3 def Av, 
Av, 
(ii) 
def 1 sin 0; — sin 6; 
Nii day ee es 2 1) ky # ke and ky,ko = 1,2,3 
(Fi, ka) sin(O;, — 9x, ) en Cos Ox, 1 ke nee 
and (iii) 


é(ky, ka) & 2) vk # ko and ky, ko = 1,2,3 
ko 


Note that: (i) Eqs.(3.54) and (3.69) > 
é, = &(2,3), é = é(3, 1) and é5 = (1,2) 


(ii) let [M(k1, k2)]’ be the transpose of M(ky, kz), then Eq.(3.68) > 


sin? (0; — 9%,) \—cos(Ox5 — 9x, ) 1 
for any (ki, ka) with ky x ko and ky, ko = 1,2,3 


TK, he) 22 [M (de, ke) |*M (Ki, ke) = —-———_—~ ( — a 


and (iii) Eqs.(3.52), (3.66), (3.67) and (3.71) > 


(el) = (Ava)? + (Avy)? = (6)'6 = [els ke)]! M(B, a) MCh, a)e(be ha) 
= [é(k1, k2)]/T' (ki, k2)]é(k1, ka), kr ke and ky, ko = 1,2,3 


where, 6¢ denotes the transpose of 6. 
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(3.61) 


(3.62) 


(3.63) 


(3.64) 


(3.65) 


(3.66) 


(3.67) 


(3.68) 


(3.69) 


(3.70) 


(3.71) 


(3.72) 


Note that: (i) Eq.(3.69) > 


(ka, ky) = P'e(ky, ka), for any (k1, ka) with ky # ko and ky, ko = 1,2,3 (3.73) 
where the 2 x 2 permutation matrix 
pee 4 (3.74) 
1 0 
has the property 
Para (3.75) 


where P~! is the inverse of P; (ii) Eq.(3.71) > 
T (ko, kt) = T(k1,k2), ker & ka and ky, kp = 1,2,3 (3.76) 
(iii) with the aid of Eqs.(3.74) and (3.75), Eq.(3.71) > 
P'T (ky, k2)P =T(k1,k2), ki # ko and hy, ka = 1,2,3 (3.77) 
and (iv) Eqs.(3.73), (3.76), and (3.77) > 


[€(ke, k1)]°T (ke, k1)]é(Ko, 1) = [€(K1, ko) |*T'(k1, k2)é(kr, ke) 


(3.78) 
ky x ko and ky, ke = 1,2,3 


Because of the relation Eq. (3.78), Eq.(3.72) represents only three independent relations. As such, with the 
aid of Eq.(3.70), one concludes that Eq.(3.72) © 


(JAVP) = (€:)'T2.3)& = (@)'TB, Ne = (6)'T(1. ee (3.79) 


By comparing Eq.(3.79) with Eq.(3.55) it is seen that the proof of Theorem 1 is completed if one can show 
that 


T(2,3) = S(a1), T(3,1) = S(a2) and T(1, 2) = S(az3) (3.80) 
To prove Eq.(3.80), note that Eqs.(3.25), (3.26), and (3.41) > 
3 PoP3. P3P, and & P,P» (3.81) 
ey = 5 62> nd €3 = i 
IPP, PsP PiP| 


In turn, with the aid of Fig.3(a) and Eq.(3.81) , one has 


— >); > 
—» + } 
€2¢€3 = eset [PP |PiPe] cosy = cosa (3.82) 
EZ © 3 ~~ To ——S]| |= = 1 : 
|PLPs| |PiP,| |PLPs| [Pi P| 


On the other hand, by using Eqs.(3.58) and elementary trigonometry, one has 


—€5 e €3 = cos(63 02) = cos(@2 03) (3.83) 
By using Eqs.(3.82) and (3.83) along with elementary trigonometry, one has 
—cos(02 — 63) = cosa, and sin? (02 — 63) = sin? ay (3.84) 


Moreover, by using Eqs.(3.60), (3.71) and (3.76), one concludes that (i) sin(@2 — 03) 4 0, and (ii) 


sin? (A5 Ao 3) 1 


T(2,3) =T(3,2) = = ( 7 1 — €08( 82 — ”) (3.85) 


Combining Eqs.(3.53), (3.84), and (3.85) one arrives at the first part of Eq.(3.80), ie., T(2,3) = S(a1). 
Moreover, because the new version of proof emerging from the above proof for T(2,3) = S(a1) 


14 of 38 


American Institute of Aeronautics and Astronautics 


after it undergoes the index permutation of Eq.(3.44) remains valid, one will arrive at the second 
part of Eq.(3.80). Similarly one can prove the third part of Eq.(3.80). Because Eqs.(3.79) and (3.80) > 
Eq.(3.55), the proof of Theorem 1 is completed. QED. 

As a preliminary for the following developments, note that, by using (i) Fig.3(a), (ii) Eq.(3.25), and (iii) 
the law of the sines, one concludes that: (i) 


sina, > 0,k = 1,2,3, for each (a1, a2,a3) € Da (3.86) 


and (ii) there exists a parameter 6 > 0 such that 


§2.3 §31 — §1,2 
sin a1 sin a2 sin a3 


where D, is defined in Eq.(3.50). Next, by using Eq.(3.40) and an equivalent of Eq.(3.87), i-e., 


= 6 > 0 for each (a1, Q2,a3) € Da (3.87) 


$23 = Bsinay, $31 = Ssinazg and s12=fsinag (6 >0 and (aj, a2,a3) € Da) (3.88) 
Eq.(3.38) can be recast as 
(sin a1 )e, + (sin ag)e€g + (sina3)e3 = 0, (a1, 2,03) € Da (3.89) 
Next, for a given (a1, Q2,a3) € Da, let 


A(ayz, @2,03) = {(€1, €2, €3)|€1, €2 and €3 be real parameters 


3.90 
satisfying Eq.(3.89)} (a1,a2,a3) € Da oy 


ie., for a given (a1,Q2,a3) € Da, A(a1,a2,a3) contains all the elements (€),€2,¢3) which satisfy 
Eq.(3.89). In turn, with the aid of Eq.(3.86), one concludes that, for any (€1, €2,€3) € A(a,,Q@2,a3) with 
(a1, Q2,a3) € Da, any one of €1,€2, and €3, can be uniquely determined in terms of the other two by using 
Eq.(3.89). In fact, for any (€1,é2,é3) defined from any given (€1,€2,€3) € A(ai1,a2,a3) using Eq.(3.54), 
Eq.(3.89) = any one of the following six relations: 


é = J(a1,02,03)é2, 2 = [J(a1,a2,03)]~'é1 (a) 
€9 = J (a2, a3, 1 )é3, €3 = [J (az,03,01)|~ 12 (b) (a1, a2, a3) E Da (3.91) 
and €3 = J(a3, Q1,02)é1, €y = [J(a3, 01, a2)|~*é3 (c) 


Here, for any (a1, @2,a3) € Do [which © (a2, 03,01) € Da © (a3, 01, 2) € Dal, (i) 


oa se a8 
J (a1, a2, 03) = ; ( es oF =) ,(Q1, 02,03) € Da (3.92) 
sina@2 \ sinag 0 
and (ii) [J(a1, a2, a3)|~+, is explicitly given by 
Hy 1 0 sin Q, 
[J (a1, a2, a3)] aor 7 : ’ (a1, @2, a3) € Da (3.93) 
sina; \—sinag —sina3 


As a preliminary for a key future development, given any (a1,Q2,a3) € Da, next we will prove the 
following relations: 


S(a1) = [(J(ar, a2, %3 


) 
S(a2) = [(J(a2, a3, a1) 
and S(a3) = [(J(a3, a1, a2) 


)~*]' S(a2)(J(a1,02,03))* (a) 
4]! S(a3)(I(a2,a3,01))-! (b) ¢ (a1,02,03)€ Da (8.94) 
)“1]' S(a1)(J(a3,01,02))~? (c) 
Proof for Eq.(3.94): By using Eqs.(3.53) and (3.93), one has 


1 


a) ae, 
sin” a1 -sin* a2 


[(J(a1, 02, a3))~4]* S(a2)(J(a1, a2, 03))7! = 


a) : ; ‘ (3.95) 
sin* ag sin @2(sin a3 — sin a1 Cos a2) 
‘ " ; nc : . : ’ (a1,02,03) € Da 
sin @2(sina3 — sina; cosa2) sin“ a; + sinas3(sinag — 2sin a1 cos a2) 
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To simplify Eq.(3.95) note that Eqs.(3.49) and (3.50) > a3 = 7 — a1 — a@2 and thus 
sin a3 = sin(m — a, — a2) = sin(a, + a2) = sina, cos a2 + sin a2 cos ay, (3.96) 
In turn, by using Eq.(3.96) and some trignometric identities, one has 
sin a2(sin a3 — sin a1 Cos a2) = (sin? a2) cos ay (3.97) 


and 
sin? a, + sina3(sina3 — 2sina; cosaz) = sin? ag (3.98) 


Next, by substituting Eqs.(3.97) and (3.98) into Eq.(3.95) and then using Eq.(3.53) again, one arrives at 
Eq.(3.94)(a). Because Eq.(3.94)(b) is the image of Eq.(3.94)(a) under the index permutation of Eq.(3.44), by 
carrying out this permutation over the proof for Eq.(3.94)(a), one has the proof for Eq.(3.94)(b). Similarly 
one can prove Eq.(3.94)(c). QED. 

Note that an immediate result of Eq.(3.55) is 


(€)'S(a1)& = (€2)'S(az)é = (é3)'S(ag)és, (a1, 2,03) € Da (3.99) 


In the following, Eq.(3.99) is proved directly using Eqs.(3.91) and (3.94). 
Proof for Eq.(3.99): By using the second part of Eq.(3.91)(a), one has 


(€)'S(ag)éo = (é,) [(J(ar, Q2, a3))~*] ‘ S(az2)(J (a1, Q2, a3)) te, (3.100) 
In turn, by substituting Eq.(3.94)(a) into Eq.(3.100), one has 
(€2)'S(az)ée = (€1)'S(ay)éy (3.101) 


ie., the validity of the first equality sign in Eq.(3.99) has been proved. Similary, one can prove the validity 
of the second equality sign. QED. 

Moreover, Eq.(3.91) implies that any (€1,€2,€3) which are defined in terms of the same (€1,€2,€3) € 
A(a1,@2,@3) using Eq.(3.54) has the property that é, = 0 for any k = 1,2,3 6 & =€&=é; =0. As 
such, any (€1, €2,€3) defined using the same (€1, €2,€3) € A(a1,@2,a3) must meet one of the following two 
mutually exclusive conditions: 


(a) €, = 0, for each k = 1,2,3 (3.102) 


and 
(b) é, #0, ie., (€,)°€, > 0, for each k = 1,2,3 (3.103) 


With the aid of Eq.(3.54), it is seen that Eq.(3.102) = 
€& = €&g=—=€3 = 0 (3.104) 


According to Eqs.(3.35)—(3.37) and (3.40), Eq.(3.104) implies that the numerical errors of the directional 
derivatives of ¢ evaluated along the three sides of AP, P2P3 are all zero. 
Because (i) Eq.(3.102) <= Eq.(3.104), (ii) Eq.(3.104) and the condition 


et S Va)? + (a)? + (63)? > 0 (3.105) 
are mutually exclusive, and (iii) any (€1, €2, €3) defined in terms of the same (€1, €2,€3) € A(a1, @2, a3) must 
belong to one of the two exclusive cases, Eqs.(3.102) and (3.103), one concludes that Eq.(3.103) <= Eq.(3.105). 

Because the case Eq.(3.104) is numerically unrealistic and, according to Eqs.(3.54) and (3.55), it leads 
to the trivial result Av, = Av, = 0 (which, according to Eqs.(3.40)—(3.43) and (3.45)—(3.47), simply states 
that the numerical error of the vector Vo on the plane [ is zero if the numerical errors €1, €2 and €3 of the 
directional derivatives of ¢ evaluated along the three sides of AP; P2 P3 are all zero), hereafter, for any given 
(a1, Q2,a3) € Da, unless specified otherwise, we will consider only the case Eq.(3.103), i.e., only the case 


(€1, €2, €3) € A'(a1, a2, 03), (a1, @2, a3) € Da (3.106) 
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where, 


A’ (a1, 02, 03) e {(€1, €2, €3)|(€1, €2, €3) € A(a1, a2, a3) and also satisfies 


3.107 
the condition Eq.(3.103) or its equivalent «* > 0} for any (a1, a2,a3) € Da ( ) 
Next, with the aid of Eqs.(3.86), (3.89) and (3.105), we will prove the following theorem: 
Theorem 2: Let (a) (a1,a2,a3) € Da, and thus Eq.(3.86) > 
Gt Ea aS (3.108) 
sin ay, 
(b) 
i= sie a ees (3.109) 
2°C13 14+ (C13)? 
pan eo ea (3.110) 
C23-Co1 14+ (C21) 
and 
+ (C31) 1 ? C31 ‘ C39 (3.111) 
1° C3, 2 1+ (C39)? 
(c) 
Mie S14 (Cra)? + (Crs)?, Ae 21 (3.112) 
doy 21+ (Cog)? + (Co1)?, Xo SZ 1 (3.113) 
moe def def 
Ase = 1+ (C31)? + (C3,2)*, Az = 1 (3.114) 
(d) 
Ww, = : VA C12 VAr+ C13 (3.115) 
J (Ci,2)? + (C13) C13 —C1,2 
Ww, = : VAz+ C23 V Aart C21 (3.116) 
V (C2,3)? + (C2,1) C; 1 —C,3 
and 
Ws def /A34 C31 A384 C32 (3.117) 
JV (C31)? + (C3,2) C39 —C31 
and (e) 
dy Weer, k= 1,2,3 (3.118) 


where €,,k4 = 1,2,3, are defined by Eq.(3.54) and by assumption, belong to the case Eq.(3.103). Then we 
have: 
(a) The real symmetric matrices E,, 2 and E3 are all positive definite [Ref.9, p.250]. 


(b) 


(e*)? = (é,)'Exex > 0 for each k = 1,2,3 and each (€1, €2,€3) € A’ (a1, a2, a3) (3.119) 
(c) 
Ey = (We)'We, k = 1,2,3 (3.120) 
and (d) 
(e*)? = (Wp)'dx > 0 for each k = 1,2,3, and each (€1, €2,€3) € A’ (a1, 2, 03) (3.121) 


Proof: Note that, for each k = 1,2,3,A,4 and A,_ defined in one of Eqs.(3.112)—(3.114) are the eigen- 
values of the real symmetric matrix E;,, defined in one of Eqs.(3.109)—(3.111). On the other hand, according 
to Eqs.(3.108) and (3.112)—(3.114), we have 


Art > Ar—- = 1, for each k = 1, 2,3 (3.122) 
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Thus, for each k = 1, 2,3, the eigenvalues of EF; are positive. As a result, by definition [Ref.9, p.250], for each 
k = 1,2,3, the real symmetric matrix EF, is positive definite. In turn, this coupled with Eq.(3.103) implies 
that [Ref.9, p.250]: 


(€,)'E,é, > 0 for each k = 1,2,3, and each (€1, €2,€3) € A’(a4, a2, 03) (3.123) 
Next, note that Eqs.(3.89) and (3.108) > 
€) = —(Ci2€2 + C13 €3), with Ci. > 0 and C3 > 0 (3.124) 
Then, by using Eqs.(3.54), (3.105), (3.109), (3.123) and (3.124), one concludes that 
(e*)? = [1 + (C12)*](€2)? + [1 + (C1,3)7] (€s)? + 2C1,2 « Cr, - (€2es) 


2 \tTP 2 ; (3.125) 
= (€,)"E1€, > 0 for each (€1, €2, €3) € A’(a1, a2, a3) 


i.e., the & = 1 case of Eq.(3.119) has been proved. To prove the k = 1 case of Eqs.(3.120) and (3.121), note 
that (i) 


e 1 Bw nde 1 
a, OMe) apg oe ae C13 (3.126) 
V(C1,2)? + (C13)? \C13 V(C1,2)? + (C13)? \-C1,2 
are the eigenvectors of E' with eigenvalues 41, and »,_ respectively; and (ii) 
(€14)°14 = (1_)*a_ = 1 (3.127) 


Moreover, as expected from Eq.(3.122), and a matrix theorem, i.e., two eigenvectors of a real symmetric 
matrix with different eigenvalues must be orthogonal to each other [Ref.9, p.222], we have 


(Gia) Pee = (6) E14: = 0 (3.128) 


Because €;, and €,_ satisfy Eqs. (3.127) and (3.128), by definition [Ref.9, p.122], they form a pair of 
orthonormal eigenvectors of the matrix F}. 


Let 
= : oe 8 (3.129) 
V(C1,2)? + (C13)? \Ci3 —C12 


i.e., the eigenvector €;4 and €;_ of FE, respectively, are the first and second column of matrix U;. By using 
orthorormal property of €14 and €;_, it can be shown that 


2 


(U,)*U, = U,(U;)* = (U1)? = Ip © the 2 x 2 identity matrix (3.130) 


Thus U, is a real symmetric orthogonal matrix [Ref.9, p.126] satisfying the relation 
(i) =) = (3.131) 


Moreover, because U; is formed using two orthonormal eigenvectors €;4 and €;_ of FE), the latter can be 
diagonalized through a similarity transformation involving U, [Ref.9, p.223], i.e., 


os A, 0 My 0 
U,) (EU, = = 3.132 
( 1) 141 ( 0 a ( 0 i ( ) 


where Ay and A,_ are the eigenvalues of E, defined in Eq.(3.112). Moreover, because 414 > A1— = 1 and 
thus /A14 > 1, Eqs.(3.131) and (3.132) > 


E, =U, & ) (U,)-' = (U,)* ¢ eo ) U; (3.133) 


1 


Moreover, because Eqs.(3.115) and (3.129) > 


Wi = & = 4 U; (3.134) 
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Eq.(3.133) now implies the & = 1 case of Eq.(3.120), ie., 
E, =(U,) P cc ) c ie ) U, = (Wi)'W, (3.135) 


Next, by combining Eqs. (3.125), (3.135) and (3.118), one has 
(e*)? = (€:)'Eyé, = (€1)'(Wi)'Wi és = (Wié1)’'Wiér = ()*th, > 0 for the case Eq.(3.103) (3.136) 


ie, the k=1 case of Eq.(3.121) has been proved. 

To prove Egqs.(3.119)—(3.121) for the k = 2 case, note that new version of proof emerging from the above 
proof for the k=1 case after its component-identifications indices 1, 2, and 3 undergo the index permutation 
of Eq.(3.44) remain valid. In particular, after this permutation, Eqs.(3.125), (3.135) and (3.136), respectively, 
become mates: 

(e*)? = (@)' Eéy > 0, Ey = (Wo)'(W2) and (€*)? = (v2)'a2 > 0 (3.137) 


ie., the k=2, case of Eqs.(3.119)—(3.121) has been proved. Similarly, the k=3 case can be proved by carrying 
out the index permutation of Eq.(3.44) over the proof for k=2 case. As such, Theorem 2 has been proved. 
QED. 

Next, the following corollary follows directly from Eqs.(3.86), (3.108) and (3.112)-(3.118). 
Corollary to Theorem 2: Let (a 1,@2,a3) € Da. Then, (a) 


Nes inde a pas (3.138) 
sin” a, sin a, 
where 
(a1, A2, a3) def cin? ay + sin? ay + sin? ag > 0, (a1, a2,03) € Da (3.139) 


is referred to as the shape factor of AP; P2P3, hereafter; (b) 


|W,| “ the determinant of W, = — ns = — VV ec 1k=1,2,3 (3.140) 


sin Ar 


Thus, for each k = 1, 2,3, (Wx)~+ exists; (c) Eq.(3.118) = 
én = (We) de, k = 1,2,3 (3.141) 
and therefore 


Eq.(3.103) =o, 4 (") , ie, (W%)*d, > 0, for each k = 1,2,3 (3.142) 


and (d) 


C 
1,2 C13 
(Wi)* = : vo (3.143) 
V (C1,2)? + (C13)? ARES =i 


C 
2,3 C24 
(W2)* = : Vise (3.144) 
J (C23)? + (C21)? as —C2,3 


and 


C. 
a C39 
(Ws)? = : veo (3.145) 
V (03,1)? + (C3,2)? i —C31 


As a preliminary for a key future development, note that, by using Eqs.(3.93), 
(3.108), (3.115)-(3.117), (3.138)—(3.140) and (3.143)—(3.145), one has 


def 


Gy = W3 (J(a2,03,01))~'(W2)* 
= 1 ( sinag:sinag  —,/ysina, (3.146) 
(sin? ay + siti? a1) (sin? ag + sin? a1) Jysina, —sinag-sina3 
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GoW; (J(os,01,02))- (Wa)? 


= 1 6 sin a3 - sina, —,/y sin a2 (3.147) 
(sin? a3 + sin” ag)(sin? a; + sin? a2) V7sin a2 — sings - sin ay 
and 
def = = 
G3 = W2 (J(a1,02,a3)) (Wi) + 
= 1 ( sina, -sinag  —,/ysinag3 (3.148) 
(ein? a1 + sin? a3) (sin? ag + sin? a3) V7sin ag — sin 1 - sin ag 
Hereafter, only cases with (a1, a2,a3) € Da will be considered. 
Next, with the aid of Eqs.(3.146)—(3.148) and the identities: 
(sin? a2 + sin? a1)(sin? a3 + sin? a1) = sin? a2 - sin? a3 + ysin? ay (3.149) 
sin“ a3 + sin“ a2)(sin“ a] + sin* a2) = sin“ a3- sin“ a, + ysin“ a2 
~ 2 ae Pang - 2 - 2 eo ee! 3.150 
and 
(sin? ay + sin? a3)(sin? ag + sin? a3) = sin? ay - sin? ag + ysin? a3 (3.151) 
One has 
(Gu) (Gy) =I, k=1,2,3 (3.152) 


Thus, for each k = 1,2,3,G, is a real orthogonal matrix [Ref.9, p.126], i-e., 
(Gx) = (Gr)', k=1,2,3 (3.153) 
Next, by using Eqs.(3.91)(a)—(3.91)(c), (3.118), (3.141) and (3.146)—(3.148) one can show that 


be = Gtr, v3 = Gil and yy = Goib3 (3.154) 


Proof for Eq.(3.154): By using Eqs.(3.118), (3.91)(a), (3.141) and (3.148) one has 


be = Wea = Wo(J(a1, a2, 03))~1é, = Wo(J (a1, a2, 03))~ (Wi) = Ashi 


ie., the first part of Eq.(3.154) has been proved. Similary, by using Eqs.(3.118), (3.91)(b), (3.141), (3.146) 
and (3.147), one can prove the second and third parts of Eq.(3.154). QED. 
Note that, by using Eq.(3.153), it can be shown that Eq.(3.154) © 


br = (G3) ~*d2 = (Gs)'bo, be = (Gi) dbs = (Gi)*v3 and v3 = (G2) hi = (G2)*v (3.155) 
Also, by using Eqs.(3.106), (3.107) and (3.142), one concludes that, the current basic assumption that 
(€1, €2, €3) E A'(a4, a2, a3) => 


de & (") , Le, (We)'de > 0 for each k = 1,2,3 and each (a1, 02,03) € Da (3.156) 


At this juncture, also note that a result of Eq.(3.121) is 
(v)'d = (2) = (bs) (3.157) 


which can also be proved directly using Eqs.(3.153) and (3.154). Recall that Eqs.(3.153) and (3.154) are 
derived using (i) the relations given in Eqs.(3.91)(a)—(3.91)(c) which are derived from Eq.(3.89), (ii) the 
definition of Eq.(3.118) and its inverse relation Eq.(3.141), and (iii) the definitions G1,G2 and G3 given in 
Egs.(3.146)—(3.148). In other words, Eq.(3.157) represents a set of identity relations which follow directly 
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from Eq.(3.89). As such, with the aid of the fact that (€1, €2,€3) € A’(a1,a2,a3) & Eq.(3.156), Eq.(3.121) 
c an be replaced by a simpler equivalent form, i-e., 


(c*)? = (1) "thy for each Wy 4 (' 


) , ie., (d1)'d > 0 (3.158) 


Next, by using Eqs.(3.55) and (3.141), one has 
ae 
73 


5 (Wel)? (Wi) del’S (an) (We) be = (b)* Hed, k = 1,2,3 (3.159) 
where 3 
H, = 5 (Wa) 7 S(ax) (We) tb = 1,2,3 and (a1,02,03) € Da (3.160) 


According to Eq.(3.53), for each k = 1,2,3,S(a,) is a real symmetric matrix. In turn, Eq.(3.160) = for 
each k = 1, 2,3, Hz is also a real symmetric matrix. As such, for each k = 1, 2,3, the eigenvalues of H;, are 
all real [Ref.9, p.222]. Moreover, as will be shown, the matrices H,,H2 and H3 are similar [Ref.9, p.232], 
i.e., they are related by the following similarity transformations: 


Ay => (G3)-! Hy G3, Ho = (Gi)-! Az Gy and Az = (Gz)! Ay Go (3.161) 


where G1,G 2 and G3 are the matrices defined in Eqs.(3.146)—(3.148), respectively. According to a matrix 
theorem [Ref.9, p.232], similar matrices such as H,, H2 and H3 have the same eigenvalues with the same 
multiplicities. In the following, first we prove Eq.(3.161). 

Proof of Eq.(3.161): By using Eqs.(3.160), (3.94)(a), (3.148) and (3.153) one has 


Ay = 5[((Wi)*]'S(a1)(Wi)* = S[(Wi) 1] [J (a1, a2, 03))*]'S(a2)(F (a1, a2, 03))~ (Wi) 


Dl w rl w 


(Ma) (e112, 48))-"]! Wa) "Way" 52) (Wa) Wal T(ana,09))"Wa* 3 sg 


= [Wo F(a1,02,03))-*(Wi)*]¢5 (Wa) 1}'S (2) (We) [Wa F(or, 02, 08))- (1) 
= (G3)'H2G3 = (G3) H2G3 


ie., the first part of Eq.(3.161) has been proved. Next, with the aid of Eqs.(3.160), (3.94)(b), (3.146) and 
(3.153), it can be shown that the new version of Eq.(3.162) after it undergoes the index permutation of 
Eq.(3.44) is also valid, i-e., the second part of Eq.(3.161) is also valid. Similarly, once can also prove the last 
part. QED. 

Obviously, Eq.(3.161) can be cast in the following equivalent form: 


H2 = G3H,(G3)~', H3 = G,H2(G,)~! and H, = G2H3(G2)~! (3.163) 
Moreover, by using Eqs.(3.153), (3.154) and (3.163), one can show that 
(v1)! Hida = (v2)! Howe = (Ws) Has (3.164) 
In turn, with the aid of Eq.(3.164), Eq.(3.159) can be replaced by a simpler form 
3,=> rs » 
5 (lV) = (v1) ity (3.165) 


To study the eigenvalues of Hy), note that, with the aid of Eqs.(3.108) and (3.138), Eq.(3.143) > 


Sie sin a2 sin a3 
= 1 Ty i 
(Wi)-t= are RR "Get ip Os) EoD is (3.166) 
ay) - 2 sin a3 sin ag 
V sin” ag + sin” a3 Vy sina, 


In turn, by substituting Eq.(3.53), (3.139) and (3.166) into Eq.(3.160), one has 


My — : s 5 , (Q1, 2, a3) € Do (3.167) 
2 \e1 fi 
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where, 
def Sin? a2 + sin? ag + 2(cos a1)(sin a2)(sin a3) 


dy = 


Q1,02,a3) € D 
(sin? a2 + sin? a3)(sin? ay + sin? a2 + sin? a3)’ (1,2, 03) ‘ 


def COS Q1 (sin a3 — sin” a2) 


» (a1,a2,03) € Da 


sin ay (sin? a2 + sin? as) (sin? a + sin? a2 + sin? a3) 


and 
def sin? ag + sin? a3 — 2(cos a1)(sin a2) (sin a3) 


fi 


Q1,02,03) € D 
sin” a; (sin? ay + sin? a3) : (a1, 42, a3) a 


By definition, an eigenvalue o of H, is a root of the equation: 
the determinant of the matrix (H; — oIz2) =0 
Substituting Eq.(3.167) into Eq.(3.171), one has 


o (dh + fijo Fla fi (e1)?] =0 


Thus, 
O = 044 OF 0 = O1_~ 


where 


def 3 


mt = 7tatfst V(di + fi)? — Aldi - fr — (e1)?]} 


Note that, because 
(dy + fi)? — 4[di - fa — (1)7] = (di — fr)? +. 4(e1)? > 0 


(3.168) 


(3.169) 


(3.170) 


(3.171) 


(3.172) 


(3.173) 


(3.174) 


(3.175) 


the expression under the radical sign in Eq.(3.174) is non-negative, Thus, for any (a1,a2,a3) € Da, the 
eigne values 0,4 and o,_ of Hy are real always, a fact consistent with the fact that the eigenvalues of the 


real symmetric matrix H, must be real always [Ref.9, p. 222]. 
Next, with the aid of Eqs.(3.139) and (3.168)—(3.170), one has (i) 


dy- fi —(e1)? = 


(sin? ag + sin® a3)? — 4(cos? a1) (sin? a2)(sin? a3) — (cos? a1) (sin? a3 — sin? a2) 


(sin? a;)(sin? ag + sin? a3)?(sin? a; + sin? ag + sin” a3) 


(sin? a2 + sin® a3)? — (cos? a;)(sin? a2 + sin? a3)? 


(sin? a1) (sin? ag + sin? a3)?(sin? ay + sin? a2 + sin? a3) 


(sin? a1) (sin? ag + sin? a3)? 


(sin? a1) (sin? ag + sin? a3)?(sin? ay + sin? a2 + sin? a3) 
1 1 
= = 1, 02,03) € D 
(sin? a; + sin? ag +sin?a3)  Y’ (01, 02, 03) ‘i 


and (ii) 


dy defi ss 71 
7(sin? a1)(sin? a2 + sin? a3) 


where. 


qa = (sin? a1) [sin? ag + sin? a3 + 2(cos a1)(sin a2)(sin a3)] 
+ (sin? ay + sin? a2 + sin? a3)[sin? ag + sin? a3 — 2(cos a1)(sin a2)(sin a3) 
= (sin? a1) (sin? ag + sin? a3) + 2(sin? a1)(cos a1) (sin a2)(sin a3)+ 
(sin? a1 + sin? ag + sin? a3)(sin? ag + sin? a3) — 2(sin? a1)(cos a1)(sin ag) (sin a3) 
— 2(sin? ag + sin? a3)(cos a1) (sin a2)(sin a3) 
= (sin? ag + sin? a3)[2 sin? a1 + sin? ag + sin? a3 — 2(cos a1) (sin a2)(sin a3)] 


To simplify Eq.(3.178), next we will prove the identity: 


sin? ag + sin? a3 — 2(cos a1) (sin a2)(sina3) = sin? a1, (a1, a2,a3) € Da 
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(3.176) 


(3.177) 


(3.178) 


(3.179) 


Proof of Eq.(3.179): By using Eq.(3.24) and Fig.3(a), the law of cosines can be expressed as 


(81,3)? + (81,2)? = 2(cos a1)S1,3 °812= (s2,3)°, (a1, Q2, 3) E Da (3.180) 


Because 83,1 = $1,3 and o > 0, Eq.(3.179) now follows directly from Eqs.(3.88) and (3.180). Eq.(3.179) can 
also be proved directly without invoking the law of cosines. Because, (a1, 2,03) € Da > a, = 7—A2— Qs, 
we have 


sin? a, = sin?(azg + a3) = (sina2 cosa3 + sin a3 cos a2)” 


= (sin? a2) (cos” a3) + (sin? a3)(cos” az) + 2(sin az) (sin a3)(cos a2) (cos a3) 
= (sin? ag)(1 — sin* a3) + (sin? a3)(1 — sin? a2) + 2(sin a2)(sin a3)(cos a2) (cos a3) 


= sin? a2 + sin? a3 + 2(sin a2)(sin a3) cos(az + a3) = sin? ag + sin? a3 — 2(cos a1) (sin a2)(sin a3) 
i.e., Eq.(3.179) has been proved. QED. With the aid of Eq.(3.179), Eq.(3.178) now implies that 
qi = 3(sin? a1) (sin? ag + sin? a3) (3.181) 


In turn, by substituting Eq.(3.181) into Eq.(3.177), one has 


3 
dit+fi=—, (a1,02,03) € Da (3.182) 


~ 


Next, with the aid of Eqs.(3.176) and (3.182), Eqs.(3.174) and (3.175) imply that, for any given 
(a1, Q2, a3) € Dao, (i) 


o14 (01, 02,03) = o4(7) = 78 + 2,/(9/4) — 7], (01, 02,03) € Da (3.183) 
and (ii) a 
2 (4 — 7) =(di — fi)? +4(e1)? > 0, (a1, 02,03) € Da (3.184) 


In turn, Eqs.(3.139) and (3.184) > 
9 
Ss er (a1, 02,03) € Da (3.185) 


Note that, it was shown earlier that H,, H2 and H3 are similar and thus they must have the 
same eigenvalues. As such, Eq.(3.183) implies that o,(y) and o_(y) must be invariant under 
the index permutation of Eq.(3.44) — A fact that leads to the expectation that they must be 
a function of parameters (such as 7) which are invariant under the permutation. 

Moreover, it is shown in Appendix B that, by using Eq.(3.139), one can prove directly that, (i) Eq.(3.185) 
represents the range of y over Dg; (ii) 


9 
Y= 7 1 = 2 = 3 > if (a1,Q2,03) € Da (3.186) 


i.e., the shape factor y reaches its maximal value $ if and only if AP; P)P3 is equilateral; and (iii) for the 
case (@1,Q2,a3) € Da, 


> 07 & min{az, a2,a3} > 07 
a Lanes O3) (3.187) 
and max{a,, a2,03} > 7 


ie., for the case (a1, 02,03) € Da, y ~ 0F > the value of the largest internal angle approaches 7~ while 
those of the other two approach 0°. 

Moreover, by using Eqs.(3.183), (3.186) and (3.167)—(3.170), one can show that, for the case (a1, a2, a3) € 
Da, (i) 


Hy, = I for the case (a1, a2,a3) € Da (3.188) 
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(ii) 


o4(5) Pana ligne ee (3.189) 


y0t 


(iii) 


is undefined at 7 = , 


dy — 3 | 9-27+6,/(9/4) —¥ ‘Oe ; 
8 (97) 7 (9/4) — ¥ 4 
and (iv) 
. 9 
is undefined at y = ri 
ay) (3.191) 
dy) 3 [9=27-6VO/)=7) ey ey 9 
8 |) JO) —7 4 
Because, 
9 9 
9-— 27+ 6/ (9/4) -y= 4 —7)+ 274+ 6/ (9/4) —-y>0,0<7< i (3.192) 
Eq.(3.190) => 
do+(7) ‘ 9 
7s <0 if0<y¥<7 (3.193) 
On the other hand, because 
9) 2y — 6 Bay = OA BV OTD N9 = 29 + 8YV OD) 
9 — 27 + 6/ (9/4) — (3.194) 
ats sorue¢<: 
(9 — 27 + 6/(9/4) — 7) 4 
Eq.(3.191) > 
do-(y) _ ,; 9 


To study the behavior of o_(y) in the limit of y + 0*, note that (i) 


lim (3 — 2,/(9/4) —y) =0 and lim y=0 (3.196) 
yt y0+F 
and (ii) 
dee A sere : dy St, ee (3.197) 
dy Cae dy = 


Then, with the aid of Eqs.(3.196) and (3.197), an application of L’Hospital rule over the function o_ (7) 
defined in Eq.(3.183) > 


1 2 1 
lim o_(7) = : lim aa Xs= (3.198) 
yor 4 ys0+ (9/4) -—% 4 3 2 
By using Eqs.(3.189)—(3.191), (3.193), (3.195) and (3.198), one concludes that: as the value of y decreases 


from ? to 0* , (i) the value of o4(7) increases monotonically from 1 to +00; and (ii) the value of o_ (7) 
decreases monotonically from 1 towards the limit value 1/2. Thus, one has (i) 


o(yy=lso(y=ley=- (3.199) 
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1 9 
3 < 9-1) <1 <94(7) f0<y¥<q (3.200) 


1 
a(y) > ie & 04(y) 3 +0 & 7 3 :0T (3.201) 
Note that Eq.(3.139) > 


y = 2 if mar{ay, a2,a3} = . and (a1,a2,a3) € Da (3.202) 
Thus, with the aid of Eq.(3.183), one concludes that 


64(9) =04(2) = : and o_(7) = 0_(2) = : 


(3.203) 
if max{ay, a2, a3} = 3 and (a1, 2, a3) € Da 


At this juncture, note that Eqs.(3.139), (3.183), (3.185)-(3.187) and (3.199)-(3.201) imply that: (i) 
) =o_(5) =o4(Z) =1 (3.204) 
(ii) 


1 
5 < o1—(a1, 2,03) = o-(y) <1 < 014(a1, a2, a3) = o+(Y) 


a eae ; (3.205) 
if (a1, @2, a3) x (=, = =) (i.e., 0<y< -) and (a1, Q2, 43) E Da 
3°33 4 
and (iii) 
1 
a1-(1, A2,a3) = 0_(y) > Com  014(a1, 02,03) = 04(y) > +0 @ 
y 2 07 & min{a1, a2,a3} > OF and max{a1,a2,a3} > 7, (3.206) 


if (a1, 02,03) € Da 
However, by using the results presented above and the fact that an n x n real symmetric matrix possesses 
a set of n linearly independent real eigenvectors associated with its exclusive real eigenvalues [Ref.8, p.306], 
one concludes that, for each (a1, a2,a3) € Da, there exist two linearly independent eigenvectors W1_(y) and 
W14(y) such that 


Aydi_(y) =o_(y)i-(y) and Ay hg. (1) = 04. (7) 14-9) 


9 (3.207) 
for each y with 0 <7 < i and each (a1, Q@2,a3) € Da 


With the above preliminaries, we are ready to tackle the central question that motivates the current 
study. First, note that, by using Eqs.(3.51) and (3.105)-(3.107) 


def /3 \Ay| = (Ave)? + (Avy)?]/2 € € (ay, 2, 
a ES > 0, (41,6263) € A’ (a1, 02,08) (3.208) 


By definition, R is the square root of the ratio of the two simple averages, i.e., the simple average of (Av,)? 
and (Av,)?, and that of (€1)?, (eg)? and (e3)?. As such R is a measure of the relative magnitudes 
of the error norms ¢*/\/3 and \Av|/V2. According to Eqs.(3.35)—(3.37), (3.40) and (3.105), «*/V/3 
is an error norm associated with the numerical errors of the directional derivatives of the 
scalar function @ evaluated ¢ along the three sides of AP, P:P3. On the other hand, according to 
Eqgs.(3.45),(3.46) and (3.51), |Av|//2 is an error norm associated with the numerical error of the 
constant gradient vector ¢ on the plane [ in the x — y — ¢ space, defined by Eqs.(3.5) and (3.6). 
As such R is a measure of how large the error norm associated with the gradient vector ¢ on 
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T is amplified from that associated with the directional derivatives of ¢ evaluated along the 
three sides of AP, P2:P3. Thus, by its definition, R represents an error amplification factor. 

Numerical experiments reveal that the value of R could (but not necessary) become very large if AP; P2 P3 
has a large aspect ratio. In the following, we will provide a mathematically rigorous explanation of this 
pitfall and also a way to avoid it even in a case in which use of a triangular grid with a large aspect ratio is 
unavoidable. 

To proceed, note that it was shown earlier that the assumption (€1, €2,€3) € A’(a1, a2, a3) = Eq.(3.156). 
Thus, by combining Egqs.(3.158), (3.165) and (3.208), one has 


— Ps Aa 
3/2)(|Av|)? _ (w1)' Hiv 7 O\ nmi ee ee 
Fe im ( = ————— for any ¥, 4 ( ie.,(a1)’ah1 > 0) (3.209) 
(e*)? (v1 )'d1 0 
In turn, with the aid of (i) Eq.(3.209), (ii) the fact that Hy is a real symmetric matrix, and (iii) Rayleigh-Ritz 
Theorem [Ref.8, p.431], one concludes that, given any (a1, a2,a3) € Da 


o_(7) < R? <04(y) for any y # (") (3.210) 
Because, Eqs.(3.199), (3.200) and (3.208) > 
1 “ 0 
5S a_(y) <1 <0 4(4) and R> 0 for any ¥ 4 (") and any (a1,Q2,a3) € Da (3.211) 
Eq.(3.210) & 
o(y) <R< Vox(y) for any oy F (") and any (a1, 02,03) € Da (3.212) 


Moreover, Eq.(3.211) and the Rayleigh-Ritz theorem also = for any (a1, a@2,a3) € Da, 


(i) R= VJo_(y) for any y, with Hyd, = o_(y)yy and yy, 4 (") [(a1, @2,03) € Dal (3.213) 


and 


(ii) R= Jo4(9) for any yy with Hy, = 04(y)i and yy F (") [(a1, a2, 03) € Dol (3.214) 


To study the implication of Eqs.(3.212)—(3.214), first consider the special case in which ay = ag = a3 = §, 
i.e., AP, P2P3 is equilateral. According to Eqs.(3.188) and (3.199), for this case, one has 
9 9 9 T 
y= ve-() = o+(3) =landAMj=[p (a) =a2 =a03 = 3) (3.215) 
As such, for this special case, (i) Eq.(3.212) > 


R=1 for any fy, 4 (") (a1, = ag = 03 = 5) (3.216) 
and (ii) 
Ay = Lah = qr for any Yi # (") (a1 = a2 =a3 = =) (3.217) 


iLe., any wy x (") is an eigenvector of H; with unit eigenvalue if AP; P2P3 is equilateral. 
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Next consider the case in which 


(a1, Q2,a3) # ( ) and (a1, @2,a3) € Da (3.218) 


71 T 
37°33 
According to Eqs.(3.185) and (3.188), Eq.(3.218) = 


9 
0<7y< i (3.219) 


Because (i) Eq.(3.200) implies that o_(y) < o4(7) for each y with 0 < y < 7, (ii) two eigenvectors of 
Hy associated with two distinct eigenvalues, respectively, must be linearly independent [Ref.8, p.268], and 
(iii) the 2 x 2 matrix H, can have at most two linearly independent eigenvectors, one concludes that, for 
each y with 0 < 7 < 3, the eigenspace of Hj associated with either of the eigenvalues o_(y) and o+(7) is 
one-dimensional. As such, for a given 7 with 0 < 7 < 3, (i) #1, which is 2 x 1 nonzero real column matrix 
by definition, is an eigenvector of H; associated with the eigenvalue o_(y) if and only if 


; . 9 
vy =r_uwy_(7), r- £0 for any y with0 <7 < rh (3.220) 


where r_ is any real number # 0, and w_(9) is the real eigenvector of H, introduced in Eq.(3.207), and 


(ii) W1 is the real eigenvector of H, associated with the eigenvalue o,(c) if and only if 


: ; 9 
v1, =74 U14(7), r+ #0 for any y with0<y7< Z (3.221) 


where r, is any real number ¥ 0, and w1(7) is the real eigenvector of H, introduced in Eq.(3.207). In turn, 
by combining Eqs.(3.212)—(3.214) with the observations given in the above items (i) and (ii), one concludes 
that, for any y withO<y7< 2 (ii)the lower bound ,/o_(7) is attained by R if and only if v1 is in 
the form specified by Eqs.(3.220) and (ii) the upper bound of ,/o+(7) is attained by R if and 
only if 7, is in the form specified by Eqs.(3.221). As such, for the case Eq.(3.218) (or equivalently the 
case Eq.(3.219)), Eq.(3.212) implies that \/o,(y) and \/a_(7), respectively, are the least upper bound and 
the greatest lower bound of R for each y with 0 < 7 < 2. 

As noted earlier, by its definition Eq.(3.208), R is a measure of how large the numerical error norm 
associated with the gradient vector ¢ on the plane [' is amplified from that associated with the directional 
derivatives of ¢ evaluated along the three sides of AP; P2P3. On the other hand, for each given y with 
0<7< 3, Eq.(3.212) states that the value of R can vary between the greatest lower bound \/o_ (7) and the 


least upper bound \/o,(y). To minimize the possible value of the numerical error norm |Av|/v2 
associated with the gradient vector of ¢ on I’, against a given error norm ¢*/,/3 associated with 
the directional derivatives of ¢ evaluated along the three sides of AP; P)P3, it is imperative to 
minimize the least upper bound \/o+(7). 

According to the discussion given preceeding Eqs.(3.199)-(3.201), within its domain 0 < y < 2, the value 
of o+(7) increases monotonically from 1 to +00 as ¥ decreases from 3 to 0+. As such \/o+(7) (i) reaches 
its minimal value 1 if and only of y = 2, and (ii) approaches +00 if and only if in the limit of y — OT. 

Moreover, according to Eq.(3.186), y = 2 if and only if AP, P:P3 is equilateral. Thus \/o+(7) reaches 
its minimal value 1 if and only if AP; P2P3 is equilateral, a result one would expect intutively. On the other 
hand, according to Eq.(3.187), y — 0* and therefore \/o,(y) approaches +o0, if and only if the values of 
the two smaller internal angles of AP, P2P3 are both much less than one, as the case depicted in Fig.2(a). 
Luckily, even though AP; P;P3 depicted there must have a high aspect ratio, by no means the values of 
two of its internal angles must both be much less than one for a triangle with high aspect ratio. In fact, 
according to a discussion given in Sec.2, AP; P2 Ps has a high aspect ratio if and only if the value of one of 
its internal angle is much less than one. As such, a triangle with the value of only one of its internal angle 
being much less than one, such as the case depicted in Fig.2(b), is also associated with a high aspect ratio, 
In fact, according to Eq.(3.202), y = 2 for any right triangle even if the value of its smallest angle is much 
less than one and thus it has a high aspect ratio. Moreover, according to Eq.(3.203), for this case, 


the least bound of R = /o+(2) = i: = 1.22( for a right triangle). (3.222) 
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In other words, the least upper bound of R can be quiet close to its minimal value 1 even if the right triangle 
AP, P,P3 has a high aspect ratio. Thus, as long as one of its internal angles is close to the right angle, the 
least upper bound of the error amplification factor R associated with AP, P:P3 is quiet close to the minimal 
value 1 even if the triangle has a large aspect ratio. 

From the above theoretical discussions, one concludes that the accuracy of gradient calculation of a 
solution variable over a triangular spatial grid is tied to the shape factor + of individual triangles from which 
the grid is built. Specifically, if the values of y associated with the component triangles are close to its 
maximal value 2 (which is exactly attained if and only if the triangle is equilateral), the upper bounds of the 
error amplification factors R associated with the component triangles will be relatively close to its minimal 
value 1. It is also shown that the value of y associated with a triangle with a large aspect ratio can be 
close to the ideal maximal value 2 as long as one of its internal angles is close to the right triangle (ie., a 
90 degree angle). As such, accuracy deterioration of gradient calculation due to the high aspect 
ratio associated with the triangular grid used can be avoided if the grid is built from triangles 
with each of them possessing the property that one of its internal angles is close to the right 
angle. As will be shown in the following sections, the above theoretical predictions have been confirmed 
numerically. 


IV. Results and Discussion 


Numerical algorithms to circumvent numerical instability and accuracy issues for high aspect ratio tri- 
angular/tetrahedral meshes are being implemented in the software framework. The in-house CFD code 
framework ez/d® based on various schemes that belong to the CESE family has been continuously under 
development since 2004. This software framework has been written using a combined object-oriented and 
generic programming paradigms in the C++ programming language. Light-weight object-oriented hierarchy 
is used in conjunction with heavy use of template classes and functions to allow compile time polymorphism. 
Different conservation laws can be plugged in with templates that represent physics. Currently, the software 
supports either triangular/tetrahedral or quadrilateral/hexahedral unstructured meshes. Both multi-thread 
(based on low-level POSIX thread) and message passing interface (MPI) paradigms are used to facilitate 
large-scale parallel computations. Each MPI process within a computational node can be executed in multi- 
thread mode to further enhance parallel performance, especially for a memory bound multi-domain layout. 
A communication map is used for data transfer among interface zones. For a large mesh, of the order of a 
billion elements, each unstructured block can be built with its own connectivity and nodes. A global commu- 
nication map is then used to join all the independent blocks in the parallel computations. This arrangement 
allows the grid generation process to always have a low memory requirement. The interfaces among blocks 
can be continuous or discontinuous. A continuous interface mesh ensures better solution accuracy for un- 
steady flow computations. Both second- and fourth-order CESE numerical schemes are implemented for 
general conservation laws including Euler and Navier-Stokes equations in the software framework. The time 
accurate local time-stepping (TALTS) scheme,? © that allows all the elements in a mesh to march at an 
approximately uniform CFL number without violating flux conservation in time, is used to enhance parallel 
performance (smaller size elements do more computations than the larger ones) and numerical accuracy. 

High aspect ratio meshes are used routinely in viscous flow computations. In this section, several practical 
test examples are discussed to demonstrate the effectiveness of the CESE method to deal with high-aspect 
ratio meshes inside the boundary layer without any degradation in accuracy or numerical stability. The 
gradient evaluation procedures that enables CESE to use high-aspect ratio meshes are detailed in a previous 
paper.“ Furthermore, all the results featured here were computed with the use of the TALTS scheme, 
implemented in ez4d, to maintain computational efficiency and accuracy at the same time. 


IV.A. Acoustic Wave Propagation along High-Aspect Ratio Boundary-Layer Mesh 


The TALTS method is designed to handle large grid size disparities in a discretized domain. One frequently 
encountered problem in flow-physics is to simulate acoustic waves propagating through a viscous mesh. 
Computing acoustic waves using such a non-uniform mesh is quite challenging. Figure (6) shows a rectangular 
domain with a clustered triangular mesh with an aspect ratio (7) around 225 near the bottom boundary. A 
constant time step computation results in significant phase error due to the large CFL number disparities. 
The wave amplitude is also incorrectly damped. In contrast, the TALTS computation maintains a rather 
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Figure 6. Acoustic waves propagating through a domain with clustered mesh, showing density contours at a time 
instant. 


uniform CFL number (around 0.8) and thus results in much improved accuracy in both amplitude and phase 
throughout the domain. The results indicate that high aspect ratio triangular mesh poses no numerical issues 
for acoustic wave propagation. Capability to resolve acoustic or other physical waves inside the boundary 
layer is essential for direct numerical simulations of transitional or turbulent boundary layers. 


IV.B. Mach 3 Isothermal Laminar Boundary Layer 


Surface heating prediction has been one of the most important topics in many supersonic or hypersonic 
computational investigations. Accurately predicting the surface heat flux in a highly cooled boundary layer 
has important implications in vehicle thermal shield and aerodynamic design. Of interest is the isothermal 
boundary layer development behind the leading-edge bow shock of a blunt body. To isolate the boundary 
layer from shocks for numerical accuracy studies, a Mach 3 flow (with a free-stream temperature of 200 
K) over a highly cooled boundary layer with Tyyauw/Twatt,adiabatic = 0.2 is numerically computed. Due to 
the presence of the boundary layer, there is a weak leading edge shock that would make the free-stream 
conditions of the boundary layer deviate slightly from the Mach 3 conditions. The computational domain 
is a 1.2m x 0.1m rectangle (the height is about 5 boundary-layer thicknesses at the exit). A structured 
quadrilateral mesh with different sizes is sliced to form triangular elements for computations. The computed 
velocity and temperature profiles (with a mesh of 128,000 triangular elements and a maximum aspect ratio 
(n) of about 1500) at the location of Re = 10° are compared with compressible similarity boundary layer 
solutions in Fig.7, along with the surface heat flux distribution along the streamwise direction. The results 
show good agreement for laminar surface heat flux predictions. 


IV.C. RANS Computations for a Mach 2 Adiabatic Boundary Layer 


Meshes with high aspect ratio elements are most frequently used for Reynolds Averaged Navier-Stokes 
(RANS) computations. To resolve the viscous sublayer, log layer, and all the way up to the boundary layer 
edge, the near wall y* value has to be kept around or smaller than one. This step size is several orders 
of magnitude smaller than that along the strewamwise direction, which implies a very large aspect ratio. 
RANS computations are carried out using both Sparlart-Allmaras and Mentor’s SST models.!° Both skin 
friction coefficients and turbulent mean streamwise profiles are compared with results from the widely used 
NASA CFL3D code. A triangular mesh with about 105,000 elements and an aspect ratio (7) of around 3000 
obtained by slicing the quad mesh downloaded from the turbulent resource website,!° was used for a Mach 
2 adiabatic boundary layer. The agreement with structured mesh solutions is quite good (see Fig.8). It 
indicates the numerical dissipation treatment for the high aspect ratio unstructured meshes is adequate in 
resolving RANS equations. 


IV.D. Mach 11 Double Cone Benchmark Problem 


The CUBRC wind tunnel tests!! for hypersonic flows over a 15° — 25° double cone is often used to vali- 
date surface heating and the capability of CFD codes in predicting laminar or turbulent shock boundary 
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Figure 7. Comparison of Navier-Stokes solutions with compressible boundary-layer solutions at Re = 10° for a Mach 
3 flate-plate flow with Tya11/Twatl,adiabatic = 0.2 : (a) streamwise velocity (non-dimensionalized by u..) (b) temperature 
(non-dimensionalized by T,) and (c) heat flux distribution along the streamwise direction. 
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Figure 8. Comparison of RANS solutions with results from CFL3D for a Mach 2 flat-plate turbulent boundary layer 
flow with adiabatic walls: (a) skin friction coefficient (b) turbulent velocity profiles at Reg = 10,000. 
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layer interaction. Laminar calculations for the Run 35 case in Ref.[11] are carried out and the computed 
Mach number contours are shown in Fig.9 below. In all of the plots shown, the lengths have been non- 
dimensionalized with the length L of the first section of the cone. A triangular mesh with about 313,000 
elements and an aspect ratio (7) of about 400 near the wall was used for the computations. The separation 
bubble around the corner and the big subsonic bubble around the second ramp are clearly captured. The 
cone surface pressure coefficient and Stanton number are compared with the experimental data in Fig.10. 
The agreement of Cp is quite good, while the current results under-predict the surface heat flux. A grid 
sensitivity study is currently underway to sort out the heat flux discrepancies. 


Figure 9. Computed Mach number contours for the Mach 11.3 flow over the CUBRC 15-25 double cone. 
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Figure 10. Comparison of computed surface pressure coffeicent and Stanton number with experimental data for the 
Mach 11.3 flow over a double cone: (a) Cp (b) Stanton number. 


V. Conclusion 


This work discusses the fundamental accuracy issues when a mesh with very large aspect ratio elements is 
used in CFD simulations, especially inside the viscous boundary layer. Sources of inaccuracy as an outcome 
of triangular shapes are identified theoretically, followed by a discussion on potential remedies. The current 
status about the application of the CESE method in steady or unsteady computations for acoustic waves 
or flow over a viscous boundary layer with large aspect ratio triangular meshes is briefly discussed, to 
demonstrate the effectiveness of the gradient evaluation procedures to deal with mesh-related geometrical 


difficulties. 
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Appendix A 


Proposition 1 presented in Section II will be proved here. To prove Eqs.(2.20) and (2.21), first note that 
Eq.(2.17) implies that each (a%1, @%2, @%3) must meet one of the following two mutually exclusive conditions: 
T 


(a) 5 > max{ag1, M2, k3} > 0 with {ky, ke, kg} = {1, 2,3} (A.1) 


and 
(b) m > max{api, OK2, 3} > ; with {k1, ko, k3} = {1, 2,3} (A.2) 


Because, in the interval 0 < a < 7/2,sina increases monotonically from 0 to 1 as a increases 
from 0 to 7/2, it is seen that, for the case Eq.(A.1), Eq.(2.20) is a direct result of Eqs.(A.1) and (2.15). In 
the following, for the case Eq.(A.2), we will prove the Eq.(2.21) is also a result of Eqs.(2.15) and (2.17). 

First we will prove that, for the case Eq.(A.2), 


An2 F MAL{AK1, U2, KZ} and ap3 A max{ag1,A%2,A%3} (case Eq.(A.2)) (A.3) 
To prove the first part of Eq.(A.3) by contradiction, assume 
T > Apg = Max{aK1, UK2, AK3} > . and case Eq.(A.2) (A.4) 


and let te 
OR. = 7 — Ke (A.5) 


Because Eq.(A.5) & apg = 7 — Ga, Eq.(A.4) and the second part of Eq.(2.17), respectively, imply that 


> a2 > 0 for case Eq.(A.4) (A.6) 
and 
AK2 = Aki + Oks (A.7) 
In turn, Eqs.(A.6) and (A.7) > 
. > api + AK3 > 0 for case Eq.(A.4) (A.8) 


On the other hand, according to the first part of Eq.(2.17), agi > 0 and az3 > 0. Thus, Eq.(A.8) implies 
that 

> api > 0 and ; > arg > 0 for case Eq.(A.4) (A.9) 
Next, because sin(a — a) = sina for any real number a, Eq.(A.5) implies that 


sin Q,p2 = sin Apa (A.10) 
In turn, with the aid of Eq.(A.10), Eq.(2.15) implies that 
1 > sinag > sin@gg > sinaz3 > 0 (A.11) 


On the other hand, because, in the interval 0 < a < 7/2, sina increases monotonically from 0 to 1 as a 
increases from 0 to 7/2, for the current case Eqs.(A.6) and (A.9), Eq.(A.11) © 


. > ax > OR > Ong >0 (case Eqs.(A.6) and (A.9)) (A.12) 


By substituting Eq.(A.7) into a result of Eq.(A.12), we have 
Qn > AK2 = Ani +ax3 >0 (case Eqs.(A.6) and (A.9)) (A.13) 


which leads to the result 0 > a,3 and thus a contradiction to the first part of Eq.(2.17). As such the 
assumption Eq.(A.4) is false and therefore the proposition given in the first part of Eq.(A.3) is 
required for consistency to Eqs.(2.15) and (2.17). 
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Because of Eq.(A.3), for the case Eq.(A.2), consistency to Eqs.(2.15) and (2.17) requires that 
T > Ap = Max{agi, OK2, UK3} > . for case Eq.(A.2) (A.14) 


Because, 
TT TT def 
T >On > 5 > 2 OM = Tn > 0 (A.15) 


2 
one can show that, for the case Eq.(A.4), Eq.(2.17) © 


Tv ud 
a = Ol = 1 — On = One + Ons > On2 > 0 


and for case Eq.(A.14) (A.16) 


T 
3 = OR = 7 — Oni = Ona + Oks > OKs >0 
To search for the extra conditions on a1, 42 and a,z3 required for consistency to Eq.(2.15), note that, 
by using the definition of @z given above, and the relation sin(a — a) = sina for any real number a, it is 
seen that Eq.(2.15) = 
1> sin@gi > sinazg > sinaggs >0 (= Eq.(2.15)) (A.17) 


On the other hand, because (i) in the interval 0 < a < $,sina increases monotonically from 0 to 1 as a 
increases from 0 to 5, and (ii) Eq.(A.16) > 5 > agi > 0,5 > agg > Oand § > ax3 > 0, for the case 


Eq.(A.16), Eq.(A.17) & 
> We = 1 — Og > Ane > aK3 >0 case Eq.(A.16) (A.18) 


As such, for the case Eq.(A.14), Eqs.(A.16) and (A.18) = Eqs.(2.15) and (2.17). On the other hand, it can 
be shown that Eqs.(A.16) and (A.18) = 


_ > 1 — Op = A. + AK > Ane > AKZ >0 case Eq.(A.14) (A.19) 
As such, for the case Eq.(A.14), it has been shown that Eqs.(2.15) and (2.17) = Eq.(A.19). Because, it 
was shown earlier that, the case Eq.(A.14) is the only subcase of Eq.(A.2) which can be consistent with 
Eqs.(2.15) and (2.17), one concludes that, for the case Eq.(A.2), Eqs.(2.15) and (2.17) = Eq.(A.19). 
Because, Eqs.(A.19) = (i) Eq.(2.21) and (ii) the second part of Eq.(2.17), it has been proved that, for the 
case Eq.(A.2), Eq.(2.21) is a result of Eqs.(2.15) and (2.17). 

It was shown earlier that, (i) as a result of Eq.(2.17), each (a%1,x%2,%3) must meet one of the two 
mutually exclusive conditions Eqs.(A.1) and (A.2); and (ii) for the case Eq.(A.1), Eq.(2.20) is a result of 
Eq.(2.15). By combining the results with the result presented in the last paragraph, one concludes that 
Eq.(2.20) and (2.21) imply either Eq.(2.20) or Eq.(2.21). In the following, Eqs.(2.22)—(2.25) will be 
proved using Eqs.(2.16) and (2.19)—(2.21). 

Note that: (i) Eqs.(2.17) and (2.20) = 


. > ag, > axe > 0 and > a3 >0 case Eqgs.(2.17) and (2.20) (A.20) 
and (ii) Eqs.(2.17) and (2.21) 
5 >> 1 —ax1 > ang > 0 and - > a%3 >0 case Eqs.(2.17) and (2.21) (A.21) 


Because (i) in the interval 0 < a < $,sina increases monotonically from 0 to 1 as a increases from 0 to 5; 
and (ii) in the interval 0 < a < 5,0 < sina < 1; and (iii) sin(w — a) = sina for any real number a, one 


concludes that (i) Eq.(A.20) > 


sin Ap] Deh V3 F 
——_ > —)=—o> AA. . 
iia 1 and sin( 3) 7 2 sinaks > 0 (case Eq.(A.20)) (A.22) 


and (ii) Eq.(A.21) > 


sina. _ sin(a — 1) 


> 1 and sin(7) = —>sinaz3>0 (case Eq.(A.21)) (A.23) 


il 
sin Ap sin ap J2 
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As such, for the case Eq.(A.20), Eqs.(2.16) and (A.22) = (i) 


1 = 2 
~ sinap3  f3 


1.155 (case Eq.(A.20)) (A.24) 


and (ii) 


2 
1 = —= if and only if ag1 = agg = ang = ; 
v3 . (A.25) 
Le., @] =Q2=a0a3 = 3 (case Eq.(A.20)) 
On the other hand, for the case Eq.(A.21), Eqs.(2.16) and (A.23) > 
1 2 
Se > V2 1.414> — (case Eq.(A.21)) (A.26) 
sin ang V3 


Eqs.(2.22) and (2.23) now follow directly from Eqs.(A.20)—(A.26). 
To prove Eqs.(2.24) and (2.25) note that both Eqs.(2.20) and (2.21) > 


5 > aK2 > an3 >0 for case Eqs.(2.20) or case Eq. (2.21) (A.27) 


By using the fact that, in the interval 0 < a < 5,cota increases monotonically from 0 to +00 as a 
decreases from 5 to 0, (i) Eq.(A.27) => 


cot ar; > cotax2 >0 for case Eqs.(2.20) or case Eq.(2.21) (A.28) 
and (ii) 
lim cota = +00 (A.29) 
aot 


On the other hand, by combining Eqs.(2.19) and (A.28), one has 
2cotar3 > > cotax3 >0 for case Eqs.(2.20) or case Eq.(2.21) (A.30) 


In turn, with the aid of Eqs.(A.27) and (A.29), Eq.(A.30) = both Eqs.(2.24) and (2.25). QED. 


Appendix B 


In this appendix, Eqs.(3.185)—(3.187) will be proved using Eq.(3.139). 
To proceed, note that, with the aid of Eq.(3.50), Eq.(3.139) > 


(a4, a2, a3) = y(a4, Q2,7 — QA, — Q2) 


def . 2 - 2 om) (B.1) 
= g(a1,Q2) = sin* a, + sin“ ag + sin“(a1 + a2), (@1,Qa2,a3) € Da 
where the domain of the function g is 
D(g) © {(a1, a2)|a1 > 0,02 > 0 and a; + a2 < 7} (B.2) 


Given Eqs.(B.1) and (B.2), next we will prove Proposition B-1: 
Proposition B-1: Within its domain D(g), the global maximum of the function g is 
at and only at 


g 


7, and it is attained 


(a1, a2) = (17/3, 7/3) (B.3) 
Proof: By using Eq.(B.1) and some trigonometric identities, one has 
Og : ; 3 
ir sin(2a1) + sin[2(a1 + a2)] = 2[sin(2a, + a2)|(cos a2) (B.4) 
a1 
Og : : ‘ 
ie sin(2a2) + sin[2(a, + a2)] = 2[sin(a, + 2a2)]|(cos a1) (B.5) 
a2 
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07g 


a 5 = 4c08(2a1 + a2)(cos a2) (B.6) 
Oa4y 
07g 
a5 = 4c0s(a1 + 2a2)(cos ay) (B.7) 
Oas5 
and 2 
Git 5 
Boas = 2 cos[2(a1 + a2)] (B.8) 


By using the last expression in each of Eq.(B.4) and (B.5), it is seen that, within D(g), the extremum 
condition [Ref.12, p. 304] 

we = SH 0, (arya) € DIV) (8.9) 
<> one of the following four possible cases: 
(a) cosa, = cosag = 0; (b) cosag = 0 and sin(a, + 2a2) = 0; (c) cosa, = 0 and sin(2a1 + ag) = 0; and 
(d) sin(2a; + az) = sin(a; + 2a2) = 0 

In D(g),0 < a1,a2 < 7. Thus case (a) & a1 = ag = 7/2 > a, + Q2 = 7. On the other hand, according 
to Eq.(B.2), a] + ag = 7 => (a1, 02) ¢ D(g). Thus case (a) cannot occur for any (a1,a2) € D(g). 

For any (a1,a2) € D(g),cosaz = 0 4 ag = 7/2. On the other hand, with ag = 5, sin(a, + 2a2) = 
sin(7 + a,) = —sina;. Because sina; > 0 for any (a1,a2) € D(g), case (b) also cannot occur for any 
(a1,Q@2) € D(g). Similarly one can show that case (c) also cannot occur for any (a@1,a2) € D(g). 

Because 2a, + a2 = (a1 + Q2) + a1 and ay + 2ag = (a1 + a2) + a2, it is seen that, for any (a1, a2) € 
D(g),0 < 2a, +a2 < 2m and 0 < a, + 2a2 < 2m. Thus case (d) © 


2a, +ag =a, +202 =7 (B.10) 


for any (a@1,@2) € D(g). Because Eq.(B.10) <= Eq.(B.3), one concludes that, within D(g), the extremum 
condition Eq.(B.9) is met if and only if a; = a2 = 3. 
Let ay = a2 = §. Then, Eqs.(B.5)-(B.7) > 
ag _ O°9 0°g 


= =-2 =-1 = = B.11 
dat Oa aad 0a, 0a2 een es9) ( ) 


In turn, Eq.(B.11) > 


0?g 07g 07g A 
— — — B.12 
daz Oa% (G24) pe AOS ree = G2) ( ) 
and 
07g 
aa = -2<0 (a1 =a2= 77/3) (B.13) 
Oa4 


when a, = a2 = 7/3. Because (i) within D(g), the extremum condition, Eq.(B.9) is met if and only if 
ay = a2 = 7/3; (ii) Eqs.(B.12) and (B.13) => a relative maximum of g occurs when a1 = a2 = 7/3 [Ref.12, 
p. 312]; (iii) D(g) is an open domain [Ref.13, p.28]; and (iv) 


9 
g(m/3, 0/3) = 7 (B.14) 
one concludes that the global maximum of g in D(g) is 9/4 and it is attained at and only at (a1,a2) = 


(x/3, 7/3). QED. 
With the aid of Proposition B-1. Eqs.(B.1) and (B.2) now imply (i) Eq.(3.186), (ii) 


0 < g(a1, a2) < 9/4 if (a1, a2) € D(g) and (a1, a2) 4 (7/3, 7/3) (B.15) 

and (iii) 
lim Q1,02) = lim 1,02) = lim a1,a2) =0 B.16 
heen : 2) eee ee 2) ae at re ‘ 2) ( ) 


Moreover, because (i) g is continuous over D(g), and (ii) points (0,0), (0,7) and (7,0) on the a; — a2 
plane (see Fig.11) are limit points [Ref.13, p.28] of D(g) one concludes from Eqs.(B.14)—(B.16) that, given 
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Ly 


1 
(0,07 “——, Nr, 0) 
Ly 


Figure 11. Representation of a triangular element in a; — az space. 


any number xz with 0 < x < 9/4, there exists an (a@2,a2) € D(g) such that g(az,a2) = x. Thus, the range 
of g over D(g) is 0 < g < 9/4. In turn, with the aid of Eqs.(B.1) and (B.2), one concludes that the range of 
y over Dg is defined by Eq.(3.185). 

To prove Eq.(3.187), first note that continuity of g over D(g) coupled with Eqs.(B.14) and (B.15) implies 
that, for any given fixed point (a?,a3) € D(g), (i) g(a?, a8) is a fixed number > 0, and (ii) 
g(a1,42) = g(a}, a3) > 0, (a1, 02) € D(g) and (a, a3) € D(g) (B.17) 


1m 
(a2,02)>(a},a9 


In turn, Eqs.(B.1), (B.2) and (B.17) imply that, as long as (a, a$) is a fixed poin € D(g), it is impossible 
that y + 0* as (a2,a2) > (a), a3) with (a1,a2) € D(g). It becomes possible only if (a9, a2) is replaced by 
some special fixed point that lies on the boundary of open domain D(g) so that it lies outside of D(g) and 
yet is a limit point of D(g). 


To search for all the special limit points referred to above, let (i) 


G(a1, 2) &' sin? ay + sin? a2 + sin?(a, +2) (a1,a2) € D(G) (B.18) 


where, 
D(G) 4 f(a1, a2)o1 > 0,a2 > 0, and a; + ag < m} (B.19) 
and (iii) the line segments £1, L2 and L3 depicted in Fig.11 be defined by 


Dy = {(a1, a2)| 0 < a1 < TT and ag = 0} (B.20) 
Lo © {(a1,02)| 0 < a2 < m and ay = 0} (B.21) 

and 
Ls © {(a1,a2)| 0 < ay < m and a1 + a2 = 7} (B.22) 


Then, (i) one can conclude from Eqs.(B.1), (B.2), (B.18) and (B.19) that the functions of g and g are identical 
except that the domain D(g) is only a proper subset of D(g); (ii) 


D(g) = D(g) U (D(9))’ (B.23) 


where (D(g))’ denotes the boundary of D(g), i.e., the set of all the limit points of the open domain D(g) 
which lie outside of it, i.e., 
D(g) N(D(g))' = 2 (B.24) 


where, @ denotes the empty set; (iii) By its definition Eq.(B.19), D(g) contains all its own limit points and 
thus it is a closed domain [Ref.13, p.28]; and (iv) with the aid of Eqs.(B.20)—(B.22), the definition of (D(g))’ 
and Fig.11, one concludes that 


(D(g))' =L1UL2ULs (B.25) 
nd 
° Ly NL = {(0,0)}, LO Ls = {(m,0)}, and LoN Ls = {(0,7)} (B.26) 
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Next, let 
s(a) “2sin2a, O<a<r (B.27) 
Then, by using Eqs.(B.18)—(B.22), one has 


g(a1, @2) = s(a1), 0O<a, <rif (Q1, Q2) Ely (B.28) 

G(a1, @2) = s(a2), 0<ag<7if (a1, Q2) € Ly (B.29) 
and 

9(a1, @2) = s(ai), O< a1 < wif (a1, a2) € Ls (B.30) 


On the other hand, Eq.(B.27) and the definitions 


def ds def ad? Ss 


s(a) = rs and s(a) = do2’ O<a<a (B.31) 
=> 
$(q@) = 2sin(2a@) and &(a) = 4cos(2a), O0<a<a (B.32) 
Eqs.(B.32), in turn, implies that 
Sa) =O0Sa=7/2 if0<a<r (B.33) 
and 
§(7/2) = -4 <0 (B.34) 


By combining Eqs.(B.27), (B.33) and (B.34), one concludes that, (i) no local (relative) minimum of s exists 
in the interval 0 < a < 7, (ii) the global maximal value of s in the interval 0 < a < 7 is s(m/2) = 2, and it 
is attained if and only if a = 7/2, ie., 


s(7/2) =2>s(a) >0if0<a<a7/20r7t/2<a<a7 (B.35) 


and (iii) 

s(0) = s(t) =0< s(a), O<a<a7 (B.36) 
With the aid of the results presented in the above items (i)—(iii), Eqs.(B.25) and (B.28)—Eqs.(B.30) imply 
that, over the boundary (D(g))’ of the open domain D(g), (i) the maximal value of g is 2, and 
it is attained at and only at (a), a2) = (7/2,0), or (a1,a2) = (0,7/2), or (a1, a2) = (7/2, 7/2); and 
(ii) the minimal value of g is 0, and it is attained at and only at (a1, a2) = (0,0) or (a1, a2) = 
(7,0) or (a1, a2) = (0,7). Because (a) the function g is continuous over its domain D(g), and (b) the points 
(0,0), (7,0) and (0,7) on the aj — a2 plane are the only limit points of the open domain D(g) which lie in 
(D(g))’ and yet the function g attains its minimal value 0 at these points, Eq.(3.187) now follows directly 
from Eqs.(B.1), (B.2), (B.15)-(B.19) and (B.23). Note that the fact that the maximal value of g attained in 
(D(g))’ is 2, is consistent with Eq.(3.186). 
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